ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.18
Committed: Tue Aug 26 19:10:10 2014 UTC (9 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad4R2P1
Changes since 2.17: +3 -3 lines
Log Message:
Fixed Reinhart sky sender sampling

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rfluxmtx.c,v 2.17 2014/08/26 18:56:19 greg Exp $";
3 #endif
4 /*
5 * Calculate flux transfer matrix or matrices using rcontrib
6 */
7
8 #include "copyright.h"
9
10 #include <ctype.h>
11 #include <stdlib.h>
12 #include "rtio.h"
13 #include "rtmath.h"
14 #include "rtprocess.h"
15 #include "bsdf.h"
16 #include "bsdf_m.h"
17 #include "random.h"
18 #include "triangulate.h"
19 #include "platform.h"
20
21 #ifdef getc_unlocked /* avoid horrendous overhead of flockfile */
22 #undef getc
23 #define getc getc_unlocked
24 #endif
25
26 #ifdef _WIN32
27 #define SPECIALS " \t\"$*?"
28 #define QUOTCHAR '"'
29 #else
30 #define SPECIALS " \t\n'\"()${}*?[];|&"
31 #define QUOTCHAR '\''
32 #define ALTQUOT '"'
33 #endif
34
35 #define MAXRCARG 512
36
37 char *progname; /* global argv[0] */
38
39 int verbose = 0; /* verbose mode (< 0 no warnings) */
40
41 char *rcarg[MAXRCARG+1] = {"rcontrib", "-fo+"};
42 int nrcargs = 2;
43
44 const char overflowerr[] = "%s: too many arguments!\n";
45
46 #define CHECKARGC(n) if (nrcargs >= MAXRCARG-(n)) \
47 { fprintf(stderr, overflowerr, progname); exit(1); }
48
49 int sampcnt = 0; /* sample count (0==unset) */
50
51 char *reinhfn = "reinhartb.cal";
52 char *shirchiufn = "disk2square.cal";
53 char *kfullfn = "klems_full.cal";
54 char *khalffn = "klems_half.cal";
55 char *kquarterfn = "klems_quarter.cal";
56
57 /* string indicating parameters */
58 const char PARAMSTART[] = "@rfluxmtx";
59
60 /* surface type IDs */
61 #define ST_NONE 0
62 #define ST_POLY 1
63 #define ST_RING 2
64 #define ST_SOURCE 3
65
66 typedef struct surf_s {
67 struct surf_s *next; /* next surface in list */
68 void *priv; /* private data (malloc'ed) */
69 char sname[32]; /* surface name */
70 FVECT snrm; /* surface normal */
71 double area; /* surface area / proj. solid angle */
72 short styp; /* surface type */
73 short nfargs; /* number of real arguments */
74 double farg[1]; /* real values (extends struct) */
75 } SURF; /* surface structure */
76
77 typedef struct {
78 FVECT uva[2]; /* tangent axes */
79 int ntris; /* number of triangles */
80 struct ptri {
81 float afrac; /* fraction of total area */
82 short vndx[3]; /* vertex indices */
83 } tri[1]; /* triangle array (extends struct) */
84 } POLYTRIS; /* triangulated polygon */
85
86 typedef struct param_s {
87 char hemis[32]; /* hemispherical sampling spec. */
88 int hsiz; /* hemisphere basis size */
89 int nsurfs; /* number of surfaces */
90 SURF *slist; /* list of surfaces */
91 FVECT vup; /* up vector (zero if unset) */
92 FVECT nrm; /* average normal direction */
93 FVECT udir, vdir; /* v-up tangent axes */
94 char *outfn; /* output file name (receiver) */
95 int (*sample_basis)(struct param_s *p, int, FILE *);
96 } PARAMS; /* sender/receiver parameters */
97
98 PARAMS curparams;
99 char curmod[128];
100 char newparams[1024];
101
102 typedef int SURFSAMP(FVECT, SURF *, double);
103
104 static SURFSAMP ssamp_bad, ssamp_poly, ssamp_ring;
105
106 SURFSAMP *orig_in_surf[4] = {
107 ssamp_bad, ssamp_poly, ssamp_ring, ssamp_bad
108 };
109
110 /* Clear parameter set */
111 static void
112 clear_params(PARAMS *p, int reset_only)
113 {
114 while (p->slist != NULL) {
115 SURF *sdel = p->slist;
116 p->slist = sdel->next;
117 if (sdel->priv != NULL)
118 free(sdel->priv);
119 free(sdel);
120 }
121 if (reset_only) {
122 p->nsurfs = 0;
123 memset(p->nrm, 0, sizeof(FVECT));
124 memset(p->vup, 0, sizeof(FVECT));
125 p->outfn = NULL;
126 return;
127 }
128 memset(p, 0, sizeof(PARAMS));
129 }
130
131 /* Get surface type from name */
132 static int
133 surf_type(const char *otype)
134 {
135 if (!strcmp(otype, "polygon"))
136 return(ST_POLY);
137 if (!strcmp(otype, "ring"))
138 return(ST_RING);
139 if (!strcmp(otype, "source"))
140 return(ST_SOURCE);
141 return(ST_NONE);
142 }
143
144 /* Add arguments to oconv command */
145 static char *
146 oconv_command(int ac, char *av[])
147 {
148 static char oconvbuf[2048] = "!oconv -f ";
149 char *cp = oconvbuf + 10;
150 char *recv = *av++;
151
152 if (ac-- <= 0)
153 return(NULL);
154 while (ac-- > 0) {
155 strcpy(cp, *av++);
156 while (*cp) cp++;
157 *cp++ = ' ';
158 if (cp >= oconvbuf+(sizeof(oconvbuf)-32)) {
159 fputs(progname, stderr);
160 fputs(": too many file arguments!\n", stderr);
161 exit(1);
162 }
163 }
164 strcpy(cp, recv); /* receiver goes last */
165 return(oconvbuf);
166 }
167
168 /* Check if any of the characters in str2 are found in str1 */
169 static int
170 matchany(const char *str1, const char *str2)
171 {
172 while (*str1) {
173 const char *cp = str2;
174 while (*cp)
175 if (*cp++ == *str1)
176 return(*str1);
177 ++str1;
178 }
179 return(0);
180 }
181
182
183 /* Convert a set of arguments into a command line for pipe() or system() */
184 static char *
185 convert_commandline(char *cmd, const int len, char *av[])
186 {
187 int match;
188 char *cp;
189
190 for (cp = cmd; *av != NULL; av++) {
191 const int n = strlen(*av);
192 if (cp+n >= cmd+(len-3)) {
193 fputs(progname, stderr);
194 return(NULL);
195 }
196 if ((match = matchany(*av, SPECIALS))) {
197 const int quote =
198 #ifdef ALTQUOT
199 (match == QUOTCHAR) ? ALTQUOT :
200 #endif
201 QUOTCHAR;
202 *cp++ = quote;
203 strcpy(cp, *av);
204 cp += n;
205 *cp++ = quote;
206 } else {
207 strcpy(cp, *av);
208 cp += n;
209 }
210 *cp++ = ' ';
211 }
212 if (cp <= cmd)
213 return(NULL);
214 *--cp = '\0';
215 return(cmd);
216 }
217
218 /* Open a pipe to/from a command given as an argument list */
219 static FILE *
220 popen_arglist(char *av[], char *mode)
221 {
222 char cmd[10240];
223
224 if (!convert_commandline(cmd, sizeof(cmd), av)) {
225 fputs(progname, stderr);
226 fputs(": command line too long in popen_arglist()\n", stderr);
227 return(NULL);
228 }
229 if (verbose > 0)
230 fprintf(stderr, "%s: opening pipe %s: %s\n",
231 progname, (*mode=='w') ? "to" : "from", cmd);
232 return(popen(cmd, mode));
233 }
234
235 #ifdef _WIN32
236 /* Execute system command (Windows version) */
237 static int
238 my_exec(char *av[])
239 {
240 char cmd[10240];
241
242 if (!convert_commandline(cmd, sizeof(cmd), av)) {
243 fputs(progname, stderr);
244 fputs(": command line too long in my_exec()\n", stderr);
245 return(1);
246 }
247 if (verbose > 0)
248 fprintf(stderr, "%s: running: %s\n", progname, cmd);
249 return(system(cmd));
250 }
251 #else
252 /* Execute system command in our stead (Unix version) */
253 static int
254 my_exec(char *av[])
255 {
256 char *compath;
257
258 if ((compath = getpath((char *)av[0], getenv("PATH"), X_OK)) == NULL) {
259 fprintf(stderr, "%s: cannot locate %s\n", progname, av[0]);
260 return(1);
261 }
262 if (verbose > 0) {
263 char cmd[4096];
264 if (!convert_commandline(cmd, sizeof(cmd), av))
265 strcpy(cmd, "COMMAND TOO LONG TO SHOW");
266 fprintf(stderr, "%s: running: %s\n", progname, cmd);
267 }
268 execv(compath, av); /* successful call never returns */
269 perror(compath);
270 return(1);
271 }
272 #endif
273
274 /* Get normalized direction vector from string specification */
275 static int
276 get_direction(FVECT dv, const char *s)
277 {
278 int sign = 1;
279 int axis = 0;
280
281 memset(dv, 0, sizeof(FVECT));
282 nextchar:
283 switch (*s) {
284 case '+':
285 ++s;
286 goto nextchar;
287 case '-':
288 sign = -sign;
289 ++s;
290 goto nextchar;
291 case 'z':
292 case 'Z':
293 ++axis;
294 case 'y':
295 case 'Y':
296 ++axis;
297 case 'x':
298 case 'X':
299 dv[axis] = sign;
300 return(!s[1] | isspace(s[1]));
301 default:
302 break;
303 }
304 #ifdef SMLFLT
305 if (sscanf(s, "%f,%f,%f", &dv[0], &dv[1], &dv[2]) != 3)
306 #else
307 if (sscanf(s, "%lf,%lf,%lf", &dv[0], &dv[1], &dv[2]) != 3)
308 #endif
309 return(0);
310 dv[0] *= (RREAL)sign;
311 return(normalize(dv) > 0);
312 }
313
314 /* Parse program parameters (directives) */
315 static int
316 parse_params(PARAMS *p, char *pargs)
317 {
318 char *cp = pargs;
319 int nparams = 0;
320 int i;
321
322 for ( ; ; ) {
323 switch (*cp++) {
324 case 'h':
325 if (*cp++ != '=')
326 break;
327 p->hsiz = 0;
328 i = 0;
329 while (*cp && !isspace(*cp)) {
330 if (isdigit(*cp))
331 p->hsiz = 10*p->hsiz + *cp - '0';
332 p->hemis[i++] = *cp++;
333 }
334 if (!i)
335 break;
336 p->hemis[i] = '\0';
337 p->hsiz += !p->hsiz;
338 ++nparams;
339 continue;
340 case 'u':
341 if (*cp++ != '=')
342 break;
343 if (!get_direction(p->vup, cp))
344 break;
345 while (*cp && !isspace(*cp++))
346 ;
347 ++nparams;
348 continue;
349 case 'o':
350 if (*cp++ != '=')
351 break;
352 i = 0;
353 while (*cp && !isspace(*cp++))
354 i++;
355 if (!i)
356 break;
357 *--cp = '\0';
358 p->outfn = savqstr(cp-i);
359 *cp++ = ' ';
360 ++nparams;
361 continue;
362 case ' ':
363 case '\t':
364 case '\r':
365 continue;
366 case '\n':
367 case '\0':
368 return(nparams);
369 default:
370 break;
371 }
372 break;
373 }
374 fprintf(stderr, "%s: bad parameter string: %s", progname, pargs);
375 exit(1);
376 return(-1); /* pro forma return */
377 }
378
379 /* Add receiver arguments (directives) corresponding to the current modifier */
380 static void
381 finish_receiver(void)
382 {
383 char sbuf[256];
384 int uniform = 0;
385 char *calfn = NULL;
386 char *params = NULL;
387 char *binv = NULL;
388 char *binf = NULL;
389 char *nbins = NULL;
390
391 if (!curmod[0]) {
392 fputs(progname, stderr);
393 fputs(": missing receiver surface!\n", stderr);
394 exit(1);
395 }
396 if (curparams.outfn != NULL) { /* add output file spec. */
397 CHECKARGC(2);
398 rcarg[nrcargs++] = "-o";
399 rcarg[nrcargs++] = curparams.outfn;
400 }
401 /* check arguments */
402 if (!curparams.hemis[0]) {
403 fputs(progname, stderr);
404 fputs(": missing hemisphere sampling type!\n", stderr);
405 exit(1);
406 }
407 if (normalize(curparams.nrm) == 0) {
408 fputs(progname, stderr);
409 fputs(": undefined normal for hemisphere sampling\n", stderr);
410 exit(1);
411 }
412 if (normalize(curparams.vup) == 0) {
413 if (fabs(curparams.nrm[2]) < .7)
414 curparams.vup[2] = 1;
415 else
416 curparams.vup[1] = 1;
417 }
418 /* determine sample type/bin */
419 if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1') {
420 binv = "0"; /* uniform sampling -- one bin */
421 uniform = 1;
422 } else if (tolower(curparams.hemis[0]) == 's' &&
423 tolower(curparams.hemis[1]) == 'c') {
424 /* assign parameters */
425 if (curparams.hsiz <= 1) {
426 fputs(progname, stderr);
427 fputs(": missing size for Shirley-Chiu sampling!\n", stderr);
428 exit(1);
429 }
430 calfn = shirchiufn; shirchiufn = NULL;
431 sprintf(sbuf, "SCdim=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g",
432 curparams.hsiz,
433 curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
434 curparams.vup[0], curparams.vup[1], curparams.vup[2]);
435 params = savqstr(sbuf);
436 binv = "scbin";
437 nbins = "SCdim*SCdim";
438 } else if ((tolower(curparams.hemis[0]) == 'r') |
439 (tolower(curparams.hemis[0]) == 't')) {
440 calfn = reinhfn; reinhfn = NULL;
441 sprintf(sbuf, "MF=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g",
442 curparams.hsiz,
443 curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
444 curparams.vup[0], curparams.vup[1], curparams.vup[2]);
445 params = savqstr(sbuf);
446 binv = "rbin";
447 nbins = "Nrbins";
448 } else if (tolower(curparams.hemis[0]) == 'k' &&
449 !curparams.hemis[1] |
450 (tolower(curparams.hemis[1]) == 'f') |
451 (curparams.hemis[1] == '1')) {
452 calfn = kfullfn; kfullfn = NULL;
453 binf = "kbin";
454 nbins = "Nkbins";
455 } else if (tolower(curparams.hemis[0]) == 'k' &&
456 (tolower(curparams.hemis[1]) == 'h') |
457 (curparams.hemis[1] == '2')) {
458 calfn = khalffn; khalffn = NULL;
459 binf = "khbin";
460 nbins = "Nkhbins";
461 } else if (tolower(curparams.hemis[0]) == 'k' &&
462 (tolower(curparams.hemis[1]) == 'q') |
463 (curparams.hemis[1] == '4')) {
464 calfn = kquarterfn; kquarterfn = NULL;
465 binf = "kqbin";
466 nbins = "Nkqbins";
467 } else {
468 fprintf(stderr, "%s: unrecognized hemisphere sampling: h=%s\n",
469 progname, curparams.hemis);
470 exit(1);
471 }
472 if (!uniform & (curparams.slist->styp == ST_SOURCE)) {
473 SURF *sp;
474 for (sp = curparams.slist; sp != NULL; sp = sp->next)
475 if (fabs(sp->area - PI) > 1e-3) {
476 fprintf(stderr, "%s: source '%s' must be 180-degrees\n",
477 progname, sp->sname);
478 exit(1);
479 }
480 }
481 if (calfn != NULL) { /* add cal file if needed */
482 CHECKARGC(2);
483 rcarg[nrcargs++] = "-f";
484 rcarg[nrcargs++] = calfn;
485 }
486 if (params != NULL) { /* parameters _after_ cal file */
487 CHECKARGC(2);
488 rcarg[nrcargs++] = "-p";
489 rcarg[nrcargs++] = params;
490 }
491 if (nbins != NULL) { /* add #bins if set */
492 CHECKARGC(2);
493 rcarg[nrcargs++] = "-bn";
494 rcarg[nrcargs++] = nbins;
495 }
496 if (binv != NULL) {
497 CHECKARGC(2); /* assign bin variable */
498 rcarg[nrcargs++] = "-b";
499 rcarg[nrcargs++] = binv;
500 } else if (binf != NULL) {
501 CHECKARGC(2); /* assign bin function */
502 rcarg[nrcargs++] = "-b";
503 sprintf(sbuf, "%s(%g,%g,%g,%g,%g,%g)", binf,
504 curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
505 curparams.vup[0], curparams.vup[1], curparams.vup[2]);
506 rcarg[nrcargs++] = savqstr(sbuf);
507 }
508 CHECKARGC(2); /* modifier argument goes last */
509 rcarg[nrcargs++] = "-m";
510 rcarg[nrcargs++] = savqstr(curmod);
511 }
512
513 /* Make randomly oriented tangent plane axes for given normal direction */
514 static void
515 make_axes(FVECT uva[2], const FVECT nrm)
516 {
517 int i;
518
519 uva[1][0] = 0.5 - frandom();
520 uva[1][1] = 0.5 - frandom();
521 uva[1][2] = 0.5 - frandom();
522 for (i = 3; i--; )
523 if ((-0.6 < nrm[i]) & (nrm[i] < 0.6))
524 break;
525 if (i < 0) {
526 fputs(progname, stderr);
527 fputs(": bad surface normal in make_axes!\n", stderr);
528 exit(1);
529 }
530 uva[1][i] = 1.0;
531 VCROSS(uva[0], uva[1], nrm);
532 normalize(uva[0]);
533 VCROSS(uva[1], nrm, uva[0]);
534 }
535
536 /* Illegal sender surfaces end up here */
537 static int
538 ssamp_bad(FVECT orig, SURF *sp, double x)
539 {
540 fprintf(stderr, "%s: illegal sender surface '%s'\n",
541 progname, sp->sname);
542 return(0);
543 }
544
545 /* Generate origin on ring surface from uniform random variable */
546 static int
547 ssamp_ring(FVECT orig, SURF *sp, double x)
548 {
549 FVECT *uva = (FVECT *)sp->priv;
550 double samp2[2];
551 double uv[2];
552 int i;
553
554 if (uva == NULL) { /* need tangent axes */
555 uva = (FVECT *)malloc(sizeof(FVECT)*2);
556 if (uva == NULL) {
557 fputs(progname, stderr);
558 fputs(": out of memory in ssamp_ring!\n", stderr);
559 return(0);
560 }
561 make_axes(uva, sp->snrm);
562 sp->priv = (void *)uva;
563 }
564 SDmultiSamp(samp2, 2, x);
565 samp2[0] = sqrt(samp2[0]*sp->area*(1./PI) + sp->farg[6]*sp->farg[6]);
566 samp2[1] *= 2.*PI;
567 uv[0] = samp2[0]*tcos(samp2[1]);
568 uv[1] = samp2[0]*tsin(samp2[1]);
569 for (i = 3; i--; )
570 orig[i] = sp->farg[i] + uv[0]*uva[0][i] + uv[1]*uva[1][i];
571 return(1);
572 }
573
574 /* Add triangle to polygon's list (call-back function) */
575 static int
576 add_triangle(const Vert2_list *tp, int a, int b, int c)
577 {
578 POLYTRIS *ptp = (POLYTRIS *)tp->p;
579 struct ptri *trip = ptp->tri + ptp->ntris++;
580
581 trip->vndx[0] = a;
582 trip->vndx[1] = b;
583 trip->vndx[2] = c;
584 return(1);
585 }
586
587 /* Generate origin on polygon surface from uniform random variable */
588 static int
589 ssamp_poly(FVECT orig, SURF *sp, double x)
590 {
591 POLYTRIS *ptp = (POLYTRIS *)sp->priv;
592 double samp2[2];
593 double *v0, *v1, *v2;
594 int i;
595
596 if (ptp == NULL) { /* need to triangulate */
597 ptp = (POLYTRIS *)malloc(sizeof(POLYTRIS) +
598 sizeof(struct ptri)*(sp->nfargs/3 - 3));
599 if (ptp == NULL)
600 goto memerr;
601 if (sp->nfargs == 3) { /* simple case */
602 ptp->ntris = 1;
603 ptp->tri[0].vndx[0] = 0;
604 ptp->tri[0].vndx[1] = 1;
605 ptp->tri[0].vndx[2] = 2;
606 ptp->tri[0].afrac = 1;
607 } else {
608 Vert2_list *v2l = polyAlloc(sp->nfargs/3);
609 if (v2l == NULL)
610 goto memerr;
611 make_axes(ptp->uva, sp->snrm);
612 for (i = v2l->nv; i--; ) {
613 v2l->v[i].mX = DOT(sp->farg+3*i, ptp->uva[0]);
614 v2l->v[i].mY = DOT(sp->farg+3*i, ptp->uva[1]);
615 }
616 ptp->ntris = 0;
617 v2l->p = (void *)ptp;
618 if (!polyTriangulate(v2l, add_triangle)) {
619 fprintf(stderr,
620 "%s: cannot triangulate polygon '%s'\n",
621 progname, sp->sname);
622 return(0);
623 }
624 for (i = ptp->ntris; i--; ) {
625 int a = ptp->tri[i].vndx[0];
626 int b = ptp->tri[i].vndx[1];
627 int c = ptp->tri[i].vndx[2];
628 ptp->tri[i].afrac =
629 (v2l->v[a].mX*v2l->v[b].mY -
630 v2l->v[b].mX*v2l->v[a].mY +
631 v2l->v[b].mX*v2l->v[c].mY -
632 v2l->v[c].mX*v2l->v[b].mY +
633 v2l->v[c].mX*v2l->v[a].mY -
634 v2l->v[a].mX*v2l->v[c].mY) /
635 (2.*sp->area);
636 }
637 polyFree(v2l);
638 }
639 sp->priv = (void *)ptp;
640 }
641 /* pick triangle by partial area */
642 for (i = 0; i < ptp->ntris && x > ptp->tri[i].afrac; i++)
643 x -= ptp->tri[i].afrac;
644 SDmultiSamp(samp2, 2, x/ptp->tri[i].afrac);
645 samp2[0] *= samp2[1] = sqrt(samp2[1]);
646 samp2[1] = 1. - samp2[1];
647 v0 = sp->farg + 3*ptp->tri[i].vndx[0];
648 v1 = sp->farg + 3*ptp->tri[i].vndx[1];
649 v2 = sp->farg + 3*ptp->tri[i].vndx[2];
650 for (i = 3; i--; )
651 orig[i] = v0[i] + samp2[0]*(v1[i] - v0[i])
652 + samp2[1]*(v2[i] - v0[i]) ;
653 return(1);
654 memerr:
655 fputs(progname, stderr);
656 fputs(": out of memory in ssamp_poly!\n", stderr);
657 return(0);
658 }
659
660 /* Compute sample origin based on projected areas of sender subsurfaces */
661 static int
662 sample_origin(PARAMS *p, FVECT orig, const FVECT rdir, double x)
663 {
664 static double *projsa;
665 static int nall;
666 double tarea = 0;
667 int i;
668 SURF *sp;
669 /* special case for lone surface */
670 if (p->nsurfs == 1) {
671 sp = p->slist;
672 if (DOT(sp->snrm, rdir) >= -FTINY) {
673 fprintf(stderr,
674 "%s: internal - sample behind sender '%s'\n",
675 progname, sp->sname);
676 return(0);
677 }
678 return((*orig_in_surf[sp->styp])(orig, sp, x));
679 }
680 if (p->nsurfs > nall) { /* (re)allocate surface area cache */
681 if (projsa) free(projsa);
682 projsa = (double *)malloc(sizeof(double)*p->nsurfs);
683 if (!projsa) return(0);
684 nall = p->nsurfs;
685 }
686 /* compute projected areas */
687 for (i = 0, sp = p->slist; sp != NULL; i++, sp = sp->next) {
688 projsa[i] = -DOT(sp->snrm, rdir) * sp->area;
689 tarea += projsa[i] *= (double)(projsa[i] > FTINY);
690 }
691 if (tarea <= FTINY) { /* wrong side of sender? */
692 fputs(progname, stderr);
693 fputs(": internal - sample behind all sender elements!\n",
694 stderr);
695 return(0);
696 }
697 tarea *= x; /* get surface from list */
698 for (i = 0, sp = p->slist; tarea > projsa[i]; sp = sp->next)
699 tarea -= projsa[i++];
700 return((*orig_in_surf[sp->styp])(orig, sp, tarea/projsa[i]));
701 }
702
703 /* Uniform sample generator */
704 static int
705 sample_uniform(PARAMS *p, int b, FILE *fp)
706 {
707 int n = sampcnt;
708 double samp3[3];
709 double duvw[3];
710 FVECT orig_dir[2];
711 int i;
712
713 if (fp == NULL) /* just requesting number of bins? */
714 return(1);
715
716 while (n--) { /* stratified hemisphere sampling */
717 SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
718 SDsquare2disk(duvw, samp3[1], samp3[2]);
719 duvw[2] = -sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
720 for (i = 3; i--; )
721 orig_dir[1][i] = duvw[0]*p->udir[i] +
722 duvw[1]*p->vdir[i] +
723 duvw[2]*p->nrm[i] ;
724 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
725 return(0);
726 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
727 return(0);
728 }
729 return(1);
730 }
731
732 /* Shirly-Chiu sample generator */
733 static int
734 sample_shirchiu(PARAMS *p, int b, FILE *fp)
735 {
736 int n = sampcnt;
737 double samp3[3];
738 double duvw[3];
739 FVECT orig_dir[2];
740 int i;
741
742 if (fp == NULL) /* just requesting number of bins? */
743 return(p->hsiz*p->hsiz);
744
745 while (n--) { /* stratified sampling */
746 SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
747 SDsquare2disk(duvw, (b/p->hsiz + samp3[1])/curparams.hsiz,
748 (b%p->hsiz + samp3[2])/curparams.hsiz);
749 duvw[2] = sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
750 for (i = 3; i--; )
751 orig_dir[1][i] = -duvw[0]*p->udir[i] -
752 duvw[1]*p->vdir[i] -
753 duvw[2]*p->nrm[i] ;
754 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
755 return(0);
756 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
757 return(0);
758 }
759 return(1);
760 }
761
762 /* Reinhart/Tregenza sample generator */
763 static int
764 sample_reinhart(PARAMS *p, int b, FILE *fp)
765 {
766 #define T_NALT 7
767 static const int tnaz[T_NALT] = {30, 30, 24, 24, 18, 12, 6};
768 const int RowMax = T_NALT*p->hsiz + 1;
769 const double RAH = (.5*PI)/(RowMax-.5);
770 #define rnaz(r) (r >= RowMax-1 ? 1 : p->hsiz*tnaz[r/p->hsiz])
771 int n = sampcnt;
772 int row, col;
773 double samp3[3];
774 double alt, azi;
775 double duvw[3];
776 FVECT orig_dir[2];
777 int i;
778
779 if (fp == NULL) { /* just requesting number of bins? */
780 n = 0;
781 for (row = RowMax; row--; ) n += rnaz(row);
782 return(n);
783 }
784 row = 0; /* identify row & column */
785 col = b;
786 while (col >= rnaz(row)) {
787 col -= rnaz(row);
788 ++row;
789 }
790 while (n--) { /* stratified sampling */
791 SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
792 alt = (row+samp3[1])*RAH;
793 azi = (2.*PI)*(col+samp3[2]-.5)/rnaz(row);
794 duvw[2] = cos(alt); /* measured from horizon */
795 duvw[0] = tsin(azi)*duvw[2];
796 duvw[1] = tcos(azi)*duvw[2];
797 duvw[2] = sqrt(1. - duvw[2]*duvw[2]);
798 for (i = 3; i--; )
799 orig_dir[1][i] = -duvw[0]*p->udir[i] -
800 duvw[1]*p->vdir[i] -
801 duvw[2]*p->nrm[i] ;
802 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
803 return(0);
804 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
805 return(0);
806 }
807 return(1);
808 #undef rnaz
809 #undef T_NALT
810 }
811
812 /* Klems sample generator */
813 static int
814 sample_klems(PARAMS *p, int b, FILE *fp)
815 {
816 static const char bname[4][20] = {
817 "LBNL/Klems Full",
818 "LBNL/Klems Half",
819 "INTERNAL ERROR",
820 "LBNL/Klems Quarter"
821 };
822 static ANGLE_BASIS *kbasis[4];
823 const int bi = p->hemis[1] - '1';
824 int n = sampcnt;
825 double samp2[2];
826 double duvw[3];
827 FVECT orig_dir[2];
828 int i;
829
830 if (!kbasis[bi]) { /* need to get basis, first */
831 for (i = 4; i--; )
832 if (!strcasecmp(abase_list[i].name, bname[bi])) {
833 kbasis[bi] = &abase_list[i];
834 break;
835 }
836 if (i < 0) {
837 fprintf(stderr, "%s: unknown hemisphere basis '%s'\n",
838 progname, bname[bi]);
839 return(0);
840 }
841 }
842 if (fp == NULL) /* just requesting number of bins? */
843 return(kbasis[bi]->nangles);
844
845 while (n--) { /* stratified sampling */
846 SDmultiSamp(samp2, 2, (n+frandom())/sampcnt);
847 if (!bi_getvec(duvw, b+samp2[1], kbasis[bi]))
848 return(0);
849 for (i = 3; i--; )
850 orig_dir[1][i] = duvw[0]*p->udir[i] +
851 duvw[1]*p->vdir[i] +
852 duvw[2]*p->nrm[i] ;
853 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp2[0]))
854 return(0);
855 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
856 return(0);
857 }
858 return(1);
859 }
860
861 /* Prepare hemisphere basis sampler that will send rays to rcontrib */
862 static int
863 prepare_sampler(void)
864 {
865 if (curparams.slist == NULL) { /* missing sample surface! */
866 fputs(progname, stderr);
867 fputs(": no sender surface!\n", stderr);
868 return(-1);
869 }
870 /* misplaced output file spec. */
871 if ((curparams.outfn != NULL) & (verbose >= 0))
872 fprintf(stderr, "%s: warning - ignoring output file in sender ('%s')\n",
873 progname, curparams.outfn);
874 /* check/set basis hemisphere */
875 if (!curparams.hemis[0]) {
876 fputs(progname, stderr);
877 fputs(": missing sender sampling type!\n", stderr);
878 return(-1);
879 }
880 if (normalize(curparams.nrm) == 0) {
881 fputs(progname, stderr);
882 fputs(": undefined normal for sender sampling\n", stderr);
883 return(-1);
884 }
885 if (normalize(curparams.vup) == 0) {
886 if (fabs(curparams.nrm[2]) < .7)
887 curparams.vup[2] = 1;
888 else
889 curparams.vup[1] = 1;
890 }
891 VCROSS(curparams.udir, curparams.vup, curparams.nrm);
892 if (normalize(curparams.udir) == 0) {
893 fputs(progname, stderr);
894 fputs(": up vector coincides with sender normal\n", stderr);
895 return(-1);
896 }
897 VCROSS(curparams.vdir, curparams.nrm, curparams.udir);
898 if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1')
899 curparams.sample_basis = sample_uniform;
900 else if (tolower(curparams.hemis[0]) == 's' &&
901 tolower(curparams.hemis[1]) == 'c')
902 curparams.sample_basis = sample_shirchiu;
903 else if ((tolower(curparams.hemis[0]) == 'r') |
904 (tolower(curparams.hemis[0]) == 't'))
905 curparams.sample_basis = sample_reinhart;
906 else if (tolower(curparams.hemis[0]) == 'k') {
907 switch (curparams.hemis[1]) {
908 case '1':
909 case '2':
910 case '4':
911 break;
912 case 'f':
913 case 'F':
914 case '\0':
915 curparams.hemis[1] = '1';
916 break;
917 case 'h':
918 case 'H':
919 curparams.hemis[1] = '2';
920 break;
921 case 'q':
922 case 'Q':
923 curparams.hemis[1] = '4';
924 break;
925 default:
926 goto unrecognized;
927 }
928 curparams.hemis[2] = '\0';
929 curparams.sample_basis = sample_klems;
930 } else
931 goto unrecognized;
932 /* return number of bins */
933 return((*curparams.sample_basis)(&curparams,0,NULL));
934 unrecognized:
935 fprintf(stderr, "%s: unrecognized sender sampling: h=%s\n",
936 progname, curparams.hemis);
937 return(-1);
938 }
939
940 /* Compute normal and area for polygon */
941 static int
942 finish_polygon(SURF *p)
943 {
944 const int nv = p->nfargs / 3;
945 FVECT e1, e2, vc;
946 int i;
947
948 memset(p->snrm, 0, sizeof(FVECT));
949 VSUB(e1, p->farg+3, p->farg);
950 for (i = 2; i < nv; i++) {
951 VSUB(e2, p->farg+3*i, p->farg);
952 VCROSS(vc, e1, e2);
953 p->snrm[0] += vc[0];
954 p->snrm[1] += vc[1];
955 p->snrm[2] += vc[2];
956 VCOPY(e1, e2);
957 }
958 p->area = normalize(p->snrm)*0.5;
959 return(p->area > FTINY);
960 }
961
962 /* Add a surface to our current parameters */
963 static void
964 add_surface(int st, const char *oname, FILE *fp)
965 {
966 SURF *snew;
967 int n;
968 /* get floating-point arguments */
969 if (!fscanf(fp, "%d", &n)) return;
970 while (n-- > 0) fscanf(fp, "%*s");
971 if (!fscanf(fp, "%d", &n)) return;
972 while (n-- > 0) fscanf(fp, "%*d");
973 if (!fscanf(fp, "%d", &n) || n <= 0) return;
974 snew = (SURF *)malloc(sizeof(SURF) + sizeof(double)*(n-1));
975 if (snew == NULL) {
976 fputs(progname, stderr);
977 fputs(": out of memory!\n", stderr);
978 exit(1);
979 }
980 strncpy(snew->sname, oname, sizeof(snew->sname)-1);
981 snew->sname[sizeof(snew->sname)-1] = '\0';
982 snew->styp = st;
983 snew->priv = NULL;
984 snew->nfargs = n;
985 for (n = 0; n < snew->nfargs; n++)
986 if (fscanf(fp, "%lf", &snew->farg[n]) != 1) {
987 fprintf(stderr, "%s: error reading arguments for '%s'\n",
988 progname, oname);
989 exit(1);
990 }
991 switch (st) {
992 case ST_RING:
993 if (snew->nfargs != 8)
994 goto badcount;
995 VCOPY(snew->snrm, snew->farg+3);
996 if (normalize(snew->snrm) == 0)
997 goto badnorm;
998 if (snew->farg[7] < snew->farg[6]) {
999 double t = snew->farg[7];
1000 snew->farg[7] = snew->farg[6];
1001 snew->farg[6] = t;
1002 }
1003 snew->area = PI*(snew->farg[7]*snew->farg[7] -
1004 snew->farg[6]*snew->farg[6]);
1005 break;
1006 case ST_POLY:
1007 if (snew->nfargs < 9 || snew->nfargs % 3)
1008 goto badcount;
1009 finish_polygon(snew);
1010 break;
1011 case ST_SOURCE:
1012 if (snew->nfargs != 4)
1013 goto badcount;
1014 for (n = 3; n--; ) /* need to reverse "normal" */
1015 snew->snrm[n] = -snew->farg[n];
1016 if (normalize(snew->snrm) == 0)
1017 goto badnorm;
1018 snew->area = sin((PI/180./2.)*snew->farg[3]);
1019 snew->area *= PI*snew->area;
1020 break;
1021 }
1022 if ((snew->area <= FTINY) & (verbose >= 0)) {
1023 fprintf(stderr, "%s: warning - zero area for surface '%s'\n",
1024 progname, oname);
1025 free(snew);
1026 return;
1027 }
1028 VSUM(curparams.nrm, curparams.nrm, snew->snrm, snew->area);
1029 snew->next = curparams.slist;
1030 curparams.slist = snew;
1031 curparams.nsurfs++;
1032 return;
1033 badcount:
1034 fprintf(stderr, "%s: bad argument count for surface element '%s'\n",
1035 progname, oname);
1036 exit(1);
1037 badnorm:
1038 fprintf(stderr, "%s: bad orientation for surface element '%s'\n",
1039 progname, oname);
1040 exit(1);
1041 }
1042
1043 /* Parse a receiver object (look for modifiers to add to rcontrib command) */
1044 static int
1045 add_recv_object(FILE *fp)
1046 {
1047 int st;
1048 char thismod[128], otype[32], oname[128];
1049 int n;
1050
1051 if (fscanf(fp, "%s %s %s", thismod, otype, oname) != 3)
1052 return(0); /* must have hit EOF! */
1053 if (!strcmp(otype, "alias")) {
1054 fscanf(fp, "%*s"); /* skip alias */
1055 return(0);
1056 }
1057 /* is it a new receiver? */
1058 if ((st = surf_type(otype)) != ST_NONE) {
1059 if (curparams.slist != NULL && (st == ST_SOURCE) ^
1060 (curparams.slist->styp == ST_SOURCE)) {
1061 fputs(progname, stderr);
1062 fputs(": cannot mix source/non-source receivers!\n", stderr);
1063 return(-1);
1064 }
1065 if (strcmp(thismod, curmod)) {
1066 if (curmod[0]) { /* output last receiver? */
1067 finish_receiver();
1068 clear_params(&curparams, 1);
1069 }
1070 parse_params(&curparams, newparams);
1071 newparams[0] = '\0';
1072 strcpy(curmod, thismod);
1073 }
1074 add_surface(st, oname, fp); /* read & store surface */
1075 return(1);
1076 }
1077 /* else skip arguments */
1078 if (!fscanf(fp, "%d", &n)) return(0);
1079 while (n-- > 0) fscanf(fp, "%*s");
1080 if (!fscanf(fp, "%d", &n)) return(0);
1081 while (n-- > 0) fscanf(fp, "%*d");
1082 if (!fscanf(fp, "%d", &n)) return(0);
1083 while (n-- > 0) fscanf(fp, "%*f");
1084 return(0);
1085 }
1086
1087 /* Parse a sender object */
1088 static int
1089 add_send_object(FILE *fp)
1090 {
1091 int st;
1092 char otype[32], oname[128];
1093 int n;
1094
1095 if (fscanf(fp, "%*s %s %s", otype, oname) != 2)
1096 return(0); /* must have hit EOF! */
1097 if (!strcmp(otype, "alias")) {
1098 fscanf(fp, "%*s"); /* skip alias */
1099 return(0);
1100 }
1101 /* is it a new surface? */
1102 if ((st = surf_type(otype)) != ST_NONE) {
1103 if (st == ST_SOURCE) {
1104 fputs(progname, stderr);
1105 fputs(": cannot use source as a sender!\n", stderr);
1106 return(-1);
1107 }
1108 parse_params(&curparams, newparams);
1109 newparams[0] = '\0';
1110 add_surface(st, oname, fp); /* read & store surface */
1111 return(0);
1112 }
1113 /* else skip arguments */
1114 if (!fscanf(fp, "%d", &n)) return(0);
1115 while (n-- > 0) fscanf(fp, "%*s");
1116 if (!fscanf(fp, "%d", &n)) return(0);
1117 while (n-- > 0) fscanf(fp, "%*d");
1118 if (!fscanf(fp, "%d", &n)) return(0);
1119 while (n-- > 0) fscanf(fp, "%*f");
1120 return(0);
1121 }
1122
1123 /* Load a Radiance scene using the given callback function for objects */
1124 static int
1125 load_scene(const char *inspec, int (*ocb)(FILE *))
1126 {
1127 int rv = 0;
1128 char inpbuf[1024];
1129 FILE *fp;
1130 int c;
1131
1132 if (*inspec == '!')
1133 fp = popen(inspec+1, "r");
1134 else
1135 fp = fopen(inspec, "r");
1136 if (fp == NULL) {
1137 fprintf(stderr, "%s: cannot load '%s'\n", progname, inspec);
1138 return(-1);
1139 }
1140 while ((c = getc(fp)) != EOF) { /* load receiver data */
1141 if (isspace(c)) /* skip leading white space */
1142 continue;
1143 if (c == '!') { /* read from a new command */
1144 inpbuf[0] = c;
1145 if (fgetline(inpbuf+1, sizeof(inpbuf)-1, fp) != NULL) {
1146 if ((c = load_scene(inpbuf, ocb)) < 0)
1147 return(c);
1148 rv += c;
1149 }
1150 continue;
1151 }
1152 if (c == '#') { /* parameters/comment */
1153 if ((c = getc(fp)) == EOF || ungetc(c,fp) == EOF)
1154 break;
1155 if (!isspace(c) && fscanf(fp, "%s", inpbuf) == 1 &&
1156 !strcmp(inpbuf, PARAMSTART)) {
1157 if (fgets(inpbuf, sizeof(inpbuf), fp) != NULL)
1158 strcat(newparams, inpbuf);
1159 continue;
1160 }
1161 while ((c = getc(fp)) != EOF && c != '\n')
1162 ; /* else skipping comment */
1163 continue;
1164 }
1165 ungetc(c, fp); /* else check object for receiver */
1166 c = (*ocb)(fp);
1167 if (c < 0)
1168 return(c);
1169 rv += c;
1170 }
1171 /* close our input stream */
1172 c = (*inspec == '!') ? pclose(fp) : fclose(fp);
1173 if (c != 0) {
1174 fprintf(stderr, "%s: error loading '%s'\n", progname, inspec);
1175 return(-1);
1176 }
1177 return(rv);
1178 }
1179
1180 /* Get command arguments and run program */
1181 int
1182 main(int argc, char *argv[])
1183 {
1184 char fmtopt[6] = "-faa"; /* default output is ASCII */
1185 char *xrs=NULL, *yrs=NULL, *ldopt=NULL;
1186 char *iropt = NULL;
1187 char *sendfn;
1188 char sampcntbuf[32], nsbinbuf[32];
1189 FILE *rcfp;
1190 int nsbins;
1191 int a, i;
1192 /* screen rcontrib options */
1193 progname = argv[0];
1194 for (a = 1; a < argc-2; a++) {
1195 int na;
1196 /* check for argument expansion */
1197 while ((na = expandarg(&argc, &argv, a)) > 0)
1198 ;
1199 if (na < 0) {
1200 fprintf(stderr, "%s: cannot expand '%s'\n",
1201 progname, argv[a]);
1202 return(1);
1203 }
1204 if (argv[a][0] != '-' || !argv[a][1])
1205 break;
1206 na = 1;
1207 switch (argv[a][1]) { /* !! Keep consistent !! */
1208 case 'v': /* verbose mode */
1209 verbose = 1;
1210 na = 0;
1211 continue;
1212 case 'f': /* special case for -fo, -ff, etc. */
1213 switch (argv[a][2]) {
1214 case '\0': /* cal file */
1215 goto userr;
1216 case 'o': /* force output */
1217 goto userr;
1218 case 'a': /* output format */
1219 case 'f':
1220 case 'd':
1221 case 'c':
1222 if (!(fmtopt[3] = argv[a][3]))
1223 fmtopt[3] = argv[a][2];
1224 fmtopt[2] = argv[a][2];
1225 na = 0;
1226 continue; /* will pass later */
1227 default:
1228 goto userr;
1229 }
1230 break;
1231 case 'x': /* x-resolution */
1232 xrs = argv[++a];
1233 na = 0;
1234 continue;
1235 case 'y': /* y-resolution */
1236 yrs = argv[++a];
1237 na = 0;
1238 continue;
1239 case 'c': /* number of samples */
1240 sampcnt = atoi(argv[++a]);
1241 if (sampcnt <= 0)
1242 goto userr;
1243 na = 0; /* we re-add this later */
1244 continue;
1245 case 'I': /* only for pass-through mode */
1246 case 'i':
1247 iropt = argv[a];
1248 na = 0;
1249 continue;
1250 case 'w': /* options without arguments */
1251 if (argv[a][2] != '+') verbose = -1;
1252 case 'V':
1253 case 'u':
1254 case 'h':
1255 case 'r':
1256 break;
1257 case 'n': /* options with 1 argument */
1258 case 's':
1259 case 'o':
1260 na = 2;
1261 break;
1262 case 'b': /* special case */
1263 if (argv[a][2] != 'v') goto userr;
1264 break;
1265 case 'l': /* special case */
1266 if (argv[a][2] == 'd') {
1267 ldopt = argv[a];
1268 na = 0;
1269 continue;
1270 }
1271 na = 2;
1272 break;
1273 case 'd': /* special case */
1274 if (argv[a][2] != 'v') na = 2;
1275 break;
1276 case 'a': /* special case */
1277 na = (argv[a][2] == 'v') ? 4 : 2;
1278 break;
1279 case 'm': /* special case */
1280 if (!argv[a][2]) goto userr;
1281 na = (argv[a][2] == 'e') | (argv[a][2] == 'a') ? 4 : 2;
1282 break;
1283 default: /* anything else is verbotten */
1284 goto userr;
1285 }
1286 if (na <= 0) continue;
1287 CHECKARGC(na); /* pass on option */
1288 rcarg[nrcargs++] = argv[a];
1289 while (--na) /* + arguments if any */
1290 rcarg[nrcargs++] = argv[++a];
1291 }
1292 if (a > argc-2)
1293 goto userr; /* check at end of options */
1294 sendfn = argv[a++]; /* assign sender & receiver inputs */
1295 if (sendfn[0] == '-') { /* user wants pass-through mode? */
1296 if (sendfn[1]) goto userr;
1297 sendfn = NULL;
1298 if (iropt) {
1299 CHECKARGC(1);
1300 rcarg[nrcargs++] = iropt;
1301 }
1302 if (xrs) {
1303 CHECKARGC(2);
1304 rcarg[nrcargs++] = "-x";
1305 rcarg[nrcargs++] = xrs;
1306 }
1307 if (yrs) {
1308 CHECKARGC(2);
1309 rcarg[nrcargs++] = "-y";
1310 rcarg[nrcargs++] = yrs;
1311 }
1312 if (ldopt) {
1313 CHECKARGC(1);
1314 rcarg[nrcargs++] = ldopt;
1315 }
1316 if (sampcnt <= 0) sampcnt = 1;
1317 } else { /* else in sampling mode */
1318 if (iropt) {
1319 fputs(progname, stderr);
1320 fputs(": -i, -I supported for pass-through only\n", stderr);
1321 return(1);
1322 }
1323 fmtopt[2] = (sizeof(RREAL)==sizeof(double)) ? 'd' : 'f';
1324 if (sampcnt <= 0) sampcnt = 10000;
1325 }
1326 sprintf(sampcntbuf, "%d", sampcnt);
1327 CHECKARGC(3); /* add our format & sample count */
1328 rcarg[nrcargs++] = fmtopt;
1329 rcarg[nrcargs++] = "-c";
1330 rcarg[nrcargs++] = sampcntbuf;
1331 /* add receiver arguments to rcontrib */
1332 if (load_scene(argv[a], add_recv_object) < 0)
1333 return(1);
1334 finish_receiver();
1335 if (sendfn == NULL) { /* pass-through mode? */
1336 CHECKARGC(1); /* add octree */
1337 rcarg[nrcargs++] = oconv_command(argc-a, argv+a);
1338 rcarg[nrcargs] = NULL;
1339 return(my_exec(rcarg)); /* rcontrib does everything */
1340 }
1341 clear_params(&curparams, 0); /* else load sender surface & params */
1342 if (load_scene(sendfn, add_send_object) < 0)
1343 return(1);
1344 if ((nsbins = prepare_sampler()) <= 0)
1345 return(1);
1346 CHECKARGC(3); /* add row count and octree */
1347 rcarg[nrcargs++] = "-y";
1348 sprintf(nsbinbuf, "%d", nsbins);
1349 rcarg[nrcargs++] = nsbinbuf;
1350 rcarg[nrcargs++] = oconv_command(argc-a, argv+a);
1351 rcarg[nrcargs] = NULL;
1352 /* open pipe to rcontrib process */
1353 if ((rcfp = popen_arglist(rcarg, "w")) == NULL)
1354 return(1);
1355 SET_FILE_BINARY(rcfp);
1356 #ifdef getc_unlocked
1357 flockfile(rcfp);
1358 #endif
1359 if (verbose) {
1360 fprintf(stderr, "%s: sampling %d directions", progname, nsbins);
1361 if (curparams.nsurfs > 1)
1362 fprintf(stderr, " (%d elements)\n", curparams.nsurfs);
1363 else
1364 fputc('\n', stderr);
1365 }
1366 for (i = 0; i < nsbins; i++) /* send rcontrib ray samples */
1367 if (!(*curparams.sample_basis)(&curparams, i, rcfp))
1368 return(1);
1369 return(pclose(rcfp) == 0); /* all finished! */
1370 userr:
1371 if (a < argc-2)
1372 fprintf(stderr, "%s: unsupported option '%s'", progname, argv[a]);
1373 fprintf(stderr, "Usage: %s [-v][rcontrib options] sender.rad receiver.rad [-i system.oct] [system.rad ..]\n",
1374 progname);
1375 return(1);
1376 }