ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.5
Committed: Tue Jul 22 23:21:56 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +7 -4 lines
Log Message:
Fixed some bugs... more to go

File Contents

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