ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.1
Committed: Sat Jul 19 01:19:49 2014 UTC (9 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Checking in incomplete and untested version of rfluxmtx

File Contents

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