ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.2
Committed: Mon Jul 21 15:59:47 2014 UTC (9 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +61 -69 lines
Log Message:
Getting close to working version of rfluxmtx

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rfluxmtx.c,v 2.1 2014/07/19 01:19:49 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 #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 real arguments */
72 double farg[1]; /* real values (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 *calfn = NULL;
376 char *binv = NULL;
377 char *binf = NULL;
378 char *nbins = NULL;
379 char sbuf[256];
380
381 if (!curmod[0]) {
382 fputs(progname, stderr);
383 fputs(": missing receiver surface!\n", stderr);
384 exit(1);
385 }
386 if (curparams.outfn != NULL) { /* add output file spec. */
387 CHECKARGC(2);
388 rcarg[nrcargs++] = "-o";
389 rcarg[nrcargs++] = curparams.outfn;
390 }
391 /* add bin specification */
392 if (!curparams.hemis[0]) {
393 fputs(progname, stderr);
394 fputs(": missing hemisphere sampling type!\n", stderr);
395 exit(1);
396 }
397 if (normalize(curparams.nrm) == 0) {
398 fputs(progname, stderr);
399 fputs(": undefined normal for hemisphere sampling\n", stderr);
400 exit(1);
401 }
402 if (normalize(curparams.vup) == 0)
403 if (fabs(curparams.nrm[2]) < .7)
404 curparams.vup[2] = 1;
405 else
406 curparams.vup[1] = 1;
407 if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1') {
408 binv = "0";
409 } else if (tolower(curparams.hemis[0]) == 's' &&
410 tolower(curparams.hemis[1]) == 'c') {
411 /* assign parameters */
412 if (curparams.hsiz <= 1) {
413 fputs(progname, stderr);
414 fputs(": missing size for Shirley-Chiu sampling!\n", stderr);
415 exit(1);
416 }
417 sprintf(sbuf, "SCdim=%d,Nx=%g,Ny=%g,Nz=%g,Ux=%g,Uy=%g,Uz=%g",
418 curparams.hsiz,
419 curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
420 curparams.vup[0], curparams.vup[1], curparams.vup[2]);
421 CHECKARGC(2);
422 rcarg[nrcargs++] = "-p";
423 rcarg[nrcargs++] = savqstr(sbuf);
424 calfn = shirchiufn; shirchiufn = NULL;
425 binv = "scbin";
426 nbins = "SCdim*SCdim";
427 } else if ((tolower(curparams.hemis[0]) == 'r') |
428 (tolower(curparams.hemis[0]) == 't')) {
429 calfn = reinhfn; reinhfn = NULL;
430 /* XXX Need to set number of divisions */
431 binf = "rhbin";
432 nbins = "Nrhbins";
433 } else if (tolower(curparams.hemis[0]) == 'k' &&
434 !curparams.hemis[1] |
435 (tolower(curparams.hemis[1]) == 'f') |
436 (curparams.hemis[1] == '1')) {
437 calfn = kfullfn; kfullfn = NULL;
438 binf = "kbin";
439 nbins = "Nkbins";
440 } else if (tolower(curparams.hemis[0]) == 'k' &&
441 (tolower(curparams.hemis[1]) == 'h') |
442 (curparams.hemis[1] == '2')) {
443 calfn = khalffn; khalffn = NULL;
444 binf = "khbin";
445 nbins = "Nkhbins";
446 } else if (tolower(curparams.hemis[0]) == 'k' &&
447 (tolower(curparams.hemis[1]) == 'q') |
448 (curparams.hemis[1] == '4')) {
449 calfn = kquarterfn; kquarterfn = NULL;
450 binf = "kqbin";
451 nbins = "Nkqbins";
452 } else {
453 fprintf(stderr, "%s: unrecognized hemisphere sampling: h=%s\n",
454 progname, curparams.hemis);
455 exit(1);
456 }
457 if (calfn != NULL) { /* add cal file if needed */
458 CHECKARGC(2);
459 rcarg[nrcargs++] = "-f";
460 rcarg[nrcargs++] = calfn;
461 }
462 if (nbins != NULL) { /* add #bins if set */
463 CHECKARGC(2);
464 rcarg[nrcargs++] = "-bn";
465 rcarg[nrcargs++] = nbins;
466 }
467 if (binfv != NULL) {
468 CHECKARGC(2); /* assign bin variable */
469 rcarg[nrcargs++] = "-b";
470 rcarg[nrcargs++] = binv;
471 } else if (binf != NULL) {
472 CHECKARGC(2); /* assign bin function */
473 sprintf(sbuf, "%s(%g,%g,%g,%g,%g,%g)", binf,
474 curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
475 curparams.vup[0], curparams.vup[1], curparams.vup[2]);
476 rcarg[nrcargs++] = "-b";
477 rcarg[nrcargs++] = savqstr(sbuf);
478 }
479 CHECKARGC(2); /* modifier argument goes last */
480 rcarg[nrcargs++] = "-m";
481 rcarg[nrcargs++] = savqstr(curmod);
482 }
483
484 /* Make randomly oriented tangent plane axes for given normal direction */
485 static void
486 make_axes(FVECT uva[2], const FVECT nrm)
487 {
488 int i;
489
490 uva[1][0] = 0.5 - frandom();
491 uva[1][1] = 0.5 - frandom();
492 uva[1][2] = 0.5 - frandom();
493 for (i = 3; i--; )
494 if ((-0.6 < nrm[i]) & (nrm[i] < 0.6))
495 break;
496 if (i < 0) {
497 fputs(progname, stderr);
498 fputs(": bad surface normal in make_axes!\n", stderr);
499 exit(1);
500 }
501 uva[1][i] = 1.0;
502 VCROSS(uva[0], uva[1], nrm);
503 normalize(uva[0]);
504 VCROSS(uva[1], nrm, uva[0]);
505 }
506
507 /* Illegal sender surfaces end up here */
508 static int
509 ssamp_bad(FVECT orig, SURF *sp, double x)
510 {
511 fprintf(stderr, "%s: illegal sender surface '%s'\n",
512 progname, sp->sname);
513 return(0);
514 }
515
516 /* Generate origin on ring surface from uniform random variable */
517 static int
518 ssamp_ring(FVECT orig, SURF *sp, double x)
519 {
520 FVECT *uva = (FVECT *)sp->priv;
521 double samp2[2];
522 double uv[2];
523 int i;
524
525 if (uva == NULL) { /* need tangent axes */
526 uva = (FVECT *)malloc(sizeof(FVECT)*2);
527 if (uva == NULL) {
528 fputs(progname, stderr);
529 fputs(": out of memory in ssamp_ring!\n", stderr);
530 return(0);
531 }
532 make_axes(uva, sp->snrm);
533 sp->priv = (void *)uva;
534 }
535 SDmultiSamp(samp2, 2, x);
536 samp2[0] = sp->farg[6] + sqrt(samp2[0]*sp->area*(1./PI));
537 samp2[1] *= 2.*PI;
538 uv[0] = samp2[0]*tcos(samp2[1]);
539 uv[1] = samp2[0]*tsin(samp2[1]);
540 for (i = 3; i--; )
541 orig[i] = sp->farg[i] + uv[0]*uva[0][i] + uv[1]*uva[1][i];
542 return(1);
543 }
544
545 /* Add triangle to polygon's list (call-back function) */
546 static int
547 add_triangle(const Vert2_list *tp, int a, int b, int c)
548 {
549 POLYTRIS *ptp = (POLYTRIS *)tp->p;
550 struct ptri *trip = ptp->tri + ptp->ntris++;
551
552 trip->vndx[0] = a;
553 trip->vndx[1] = b;
554 trip->vndx[2] = c;
555 return(1);
556 }
557
558 /* Generate origin on polygon surface from uniform random variable */
559 static int
560 ssamp_poly(FVECT orig, SURF *sp, double x)
561 {
562 POLYTRIS *ptp = (POLYTRIS *)sp->priv;
563 double samp2[2];
564 double *v0, *v1, *v2;
565 int i;
566
567 if (ptp == NULL) { /* need to triangulate */
568 ptp = (POLYTRIS *)malloc(sizeof(POLYTRIS) +
569 sizeof(struct ptri)*(sp->nfargs/3 - 3));
570 if (ptp == NULL)
571 goto memerr;
572 if (sp->nfargs == 3) { /* simple case */
573 ptp->ntris = 1;
574 ptp->tri[0].vndx[0] = 0;
575 ptp->tri[0].vndx[1] = 1;
576 ptp->tri[0].vndx[2] = 2;
577 ptp->tri[0].afrac = 1;
578 } else {
579 Vert2_list *v2l = polyAlloc(sp->nfargs/3);
580 if (v2l == NULL)
581 goto memerr;
582 make_axes(ptp->uva, sp->snrm);
583 for (i = v2l->nv; i--; ) {
584 v2l->v[i].mX = DOT(sp->farg+3*i, ptp->uva[0]);
585 v2l->v[i].mY = DOT(sp->farg+3*i, ptp->uva[1]);
586 }
587 ptp->ntris = 0;
588 v2l->p = (void *)ptp;
589 if (!polyTriangulate(v2l, add_triangle))
590 return(0);
591 for (i = ptp->ntris; i--; ) {
592 int a = ptp->tri[i].vndx[0];
593 int b = ptp->tri[i].vndx[1];
594 int c = ptp->tri[i].vndx[2];
595 ptp->tri[i].afrac =
596 (v2l->v[a].mX*v2l->v[b].mY -
597 v2l->v[b].mX*v2l->v[a].mY +
598 v2l->v[b].mX*v2l->v[c].mY -
599 v2l->v[c].mX*v2l->v[b].mY +
600 v2l->v[c].mX*v2l->v[a].mY -
601 v2l->v[a].mX*v2l->v[c].mY) /
602 (2.*sp->area);
603 }
604 polyFree(v2l);
605 }
606 sp->priv = (void *)ptp;
607 }
608 /* pick triangle by partial area */
609 for (i = 0; i < ptp->ntris && x > ptp->tri[i].afrac; i++)
610 x -= ptp->tri[i].afrac;
611 SDmultiSamp(samp2, 2, x/ptp->tri[i].afrac);
612 samp2[0] *= samp2[1] = sqrt(samp2[1]);
613 samp2[1] = 1. - samp2[1];
614 v0 = sp->farg + 3*ptp->tri[i].vndx[0];
615 v1 = sp->farg + 3*ptp->tri[i].vndx[1];
616 v2 = sp->farg + 3*ptp->tri[i].vndx[2];
617 for (i = 3; i--; )
618 orig[i] = v0[i] + samp2[0]*(v1[i] - v0[i])
619 + samp2[1]*(v2[i] - v0[i]) ;
620 return(1);
621 memerr:
622 fputs(progname, stderr);
623 fputs(": out of memory in ssamp_poly!\n", stderr);
624 return(0);
625 }
626
627 /* Compute sample origin based on projected areas of sender subsurfaces */
628 static int
629 sample_origin(PARAMS *p, FVECT orig, const FVECT rdir, double x)
630 {
631 static double *projsa;
632 static int nall;
633 double tarea = 0;
634 int i;
635 SURF *sp;
636 /* special case for lone surface */
637 if (p->nsurfs == 1) {
638 sp = p->slist;
639 if (DOT(sp->snrm, rdir) >= -FTINY)
640 return(0); /* behind surface! */
641 return((*orig_in_surf[sp->styp])(orig, sp, x));
642 }
643 if (p->nsurfs > nall) { /* (re)allocate surface area cache */
644 if (projsa) free(projsa);
645 projsa = (double *)malloc(sizeof(double)*p->nsurfs);
646 if (!projsa) return(0);
647 nall = p->nsurfs;
648 }
649 /* compute projected areas */
650 for (i = 0, sp = p->slist; sp != NULL; i++, sp = sp->next) {
651 projsa[i] = -DOT(sp->snrm, rdir) * sp->area;
652 tarea += projsa[i] *= (double)(projsa[i] > FTINY);
653 }
654 if (tarea <= FTINY) /* wrong side of sender? */
655 return(0);
656 tarea *= x; /* get surface from list */
657 for (i = 0, sp = p->slist; tarea > projsa[i]; sp = sp->next)
658 tarea -= projsa[i++];
659 return((*orig_in_surf[sp->styp])(orig, sp, tarea/projsa[i]));
660 }
661
662 /* Uniform sample generator */
663 static int
664 sample_uniform(PARAMS *p, int b, FILE *fp)
665 {
666 int n = sampcnt;
667 double samp3[3];
668 double duvw[3];
669 FVECT orig_dir[2];
670 int i;
671
672 if (fp == NULL) /* just requesting number of bins? */
673 return(1);
674
675 while (n--) { /* stratified hemisphere sampling */
676 SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
677 SDsquare2disk(duvw, samp3[1], samp3[2]);
678 duvw[2] = -sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
679 for (i = 3; i--; )
680 orig_dir[1][i] = duvw[0]*p->udir[i] +
681 duvw[1]*p->vdir[i] +
682 duvw[2]*p->nrm[i] ;
683 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
684 return(0);
685 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
686 return(0);
687 }
688 return(1);
689 }
690
691 /* Shirly-Chiu sample generator */
692 static int
693 sample_shirchiu(PARAMS *p, int b, FILE *fp)
694 {
695 int n = sampcnt;
696 double samp3[3];
697 double duvw[3];
698 FVECT orig_dir[2];
699 int i;
700
701 if (fp == NULL) /* just requesting number of bins? */
702 return(p->hsiz*p->hsiz);
703
704 while (n--) { /* stratified sampling */
705 SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
706 SDsquare2disk(duvw, (b/p->hsiz + samp3[1])/curparams.hsiz,
707 (b%p->hsiz + samp3[2])/curparams.hsiz);
708 duvw[2] = sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
709 for (i = 3; i--; )
710 orig_dir[1][i] = -duvw[0]*p->udir[i] -
711 duvw[1]*p->vdir[i] -
712 duvw[2]*p->nrm[i] ;
713 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
714 return(0);
715 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
716 return(0);
717 }
718 return(1);
719 }
720
721 /* Reinhart/Tregenza sample generator */
722 static int
723 sample_reinhart(PARAMS *p, int b, FILE *fp)
724 {
725 #define T_NALT 7
726 static const int tnaz[T_NALT] = {30, 30, 24, 24, 18, 12, 6};
727 const int RowMax = T_NALT*p->hsiz + 1;
728 const double RAH = (.25*PI)/(RowMax-.5);
729 #define rnaz(r) (r >= RowMax-1 ? 1 : p->hsiz*tnaz[r/p->hsiz])
730 int n = sampcnt;
731 int row, col;
732 double samp3[3];
733 double alt, azi;
734 double duvw[3];
735 FVECT orig_dir[2];
736 int i;
737
738 if (fp == NULL) { /* just requesting number of bins? */
739 n = 0;
740 for (row = RowMax; row--; ) n += rnaz(row);
741 return(n);
742 }
743 row = 0; /* identify row & column */
744 col = b;
745 while (col >= rnaz(row)) {
746 col -= rnaz(row);
747 ++row;
748 }
749 while (n--) { /* stratified sampling */
750 SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
751 alt = (row+samp3[1])*RAH;
752 azi = (2.*PI)*(col+samp3[2]-.5)/rnaz(row);
753 duvw[2] = tcos(alt); /* measured from horizon */
754 duvw[0] = tcos(azi)*duvw[2];
755 duvw[1] = tsin(azi)*duvw[2];
756 duvw[2] = sqrt(1. - duvw[2]*duvw[2]);
757 for (i = 3; i--; )
758 orig_dir[1][i] = -duvw[0]*p->udir[i] -
759 duvw[1]*p->vdir[i] -
760 duvw[2]*p->nrm[i] ;
761 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
762 return(0);
763 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
764 return(0);
765 }
766 return(1);
767 #undef rnaz
768 #undef T_NALT
769 }
770
771 /* Klems sample generator */
772 static int
773 sample_klems(PARAMS *p, int b, FILE *fp)
774 {
775 static const char bname[4][20] = {
776 "LBNL/Klems Full",
777 "LBNL/Klems Half",
778 "INTERNAL ERROR",
779 "LBNL/Klems Quarter"
780 };
781 static ANGLE_BASIS *kbasis[4];
782 const int bi = p->hemis[1] - '1';
783 int n = sampcnt;
784 double samp2[2];
785 double duvw[3];
786 FVECT orig_dir[2];
787 int i;
788
789 if (!kbasis[bi]) { /* need to get basis, first */
790 for (i = 4; i--; )
791 if (!strcasecmp(abase_list[i].name, bname[bi])) {
792 kbasis[bi] = &abase_list[i];
793 break;
794 }
795 if (i < 0) {
796 fprintf(stderr, "%s: unknown hemisphere basis '%s'\n",
797 progname, bname[bi]);
798 return(0);
799 }
800 }
801 if (fp == NULL) /* just requesting number of bins? */
802 return(kbasis[bi]->nangles);
803
804 while (n--) { /* stratified sampling */
805 SDmultiSamp(samp2, 2, (n+frandom())/sampcnt);
806 if (!bo_getvec(duvw, b+samp2[1], kbasis[bi]))
807 return(0);
808 for (i = 3; i--; )
809 orig_dir[1][i] = duvw[0]*p->udir[i] +
810 duvw[1]*p->vdir[i] +
811 duvw[2]*p->nrm[i] ;
812 if (!sample_origin(p, orig_dir[0], orig_dir[1], samp2[0]))
813 return(0);
814 if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
815 return(0);
816 }
817 return(1);
818 }
819
820 /* Prepare hemisphere basis sampler that will send rays to rcontrib */
821 static int
822 prepare_sampler(void)
823 {
824 if (curparams.slist == NULL) { /* missing sample surface! */
825 fputs(progname, stderr);
826 fputs(": no sender surface!\n", stderr);
827 return(-1);
828 }
829 if (curparams.outfn != NULL) /* misplaced output file spec. */
830 fprintf(stderr, "%s: warning - ignoring output file in sender ('%s')\n",
831 progname, curparams.outfn);
832 /* check/set basis hemisphere */
833 if (!curparams.hemis[0]) {
834 fputs(progname, stderr);
835 fputs(": missing sender sampling type!\n", stderr);
836 return(-1);
837 }
838 if (normalize(curparams.nrm) == 0) {
839 fputs(progname, stderr);
840 fputs(": undefined normal for sender sampling\n", stderr);
841 return(-1);
842 }
843 if (normalize(curparams.vup) == 0)
844 if (fabs(curparams.nrm[2]) < .7)
845 curparams.vup[2] = 1;
846 else
847 curparams.vup[1] = 1;
848 VCROSS(curparams.udir, curparams.vup, curparams.nrm);
849 if (normalize(curparams.udir) == 0) {
850 fputs(progname, stderr);
851 fputs(": up vector coincides with sender normal\n", stderr);
852 return(-1);
853 }
854 VCROSS(curparams.vdir, curparams.nrm, curparams.udir);
855 if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1')
856 curparams.sample_basis = sample_uniform;
857 else if (tolower(curparams.hemis[0]) == 's' &&
858 tolower(curparams.hemis[1]) == 'c')
859 curparams.sample_basis = sample_shirchiu;
860 else if ((tolower(curparams.hemis[0]) == 'r') |
861 (tolower(curparams.hemis[0]) == 't'))
862 curparams.sample_basis = sample_reinhart;
863 else if (tolower(curparams.hemis[0]) == 'k') {
864 switch (curparams.hemis[1]) {
865 case '1':
866 case '2':
867 case '4':
868 break;
869 case 'f':
870 case 'F':
871 case '\0':
872 curparams.hemis[1] = '1';
873 break;
874 case 'h':
875 case 'H':
876 curparams.hemis[1] = '2';
877 break;
878 case 'q':
879 case 'Q':
880 curparams.hemis[1] = '4';
881 break;
882 default:
883 goto unrecognized;
884 }
885 curparams.hemis[2] = '\0';
886 curparams.sample_basis = sample_klems;
887 } else
888 goto unrecognized;
889 /* return number of bins */
890 return((*curparams.sample_basis)(&curparams,0,NULL));
891 unrecognized:
892 fprintf(stderr, "%s: unrecognized sender sampling: h=%s\n",
893 progname, curparams.hemis);
894 return(-1);
895 }
896
897 /* Compute normal and area for polygon */
898 static int
899 finish_polygon(SURF *p)
900 {
901 const int nv = p->nfargs / 3;
902 FVECT e1, e2, vc;
903 int i;
904
905 memset(p->snrm, 0, sizeof(FVECT));
906 VSUB(e1, p->farg+3, p->farg);
907 for (i = 2; i < nv; i++) {
908 VSUB(e2, p->farg+3*i, p->farg);
909 VCROSS(vc, e1, e2);
910 p->snrm[0] += vc[0];
911 p->snrm[1] += vc[1];
912 p->snrm[2] += vc[2];
913 VCOPY(e1, e2);
914 }
915 p->area = normalize(p->snrm)*0.5;
916 return(p->area > FTINY);
917 }
918
919 /* Add a surface to our current parameters */
920 static void
921 add_surface(int st, const char *oname, FILE *fp)
922 {
923 SURF *snew;
924 int n;
925 /* get floating-point arguments */
926 if (!fscanf(fp, "%d", &n)) return;
927 while (n-- > 0) fscanf(fp, "%*s");
928 if (!fscanf(fp, "%d", &n)) return;
929 while (n-- > 0) fscanf(fp, "%*d");
930 if (!fscanf(fp, "%d", &n) || n <= 0) return;
931 snew = (SURF *)malloc(sizeof(SURF) + sizeof(double)*(n-1));
932 if (snew == NULL) {
933 fputs(progname, stderr);
934 fputs(": out of memory!\n", stderr);
935 exit(1);
936 }
937 strncpy(snew->sname, oname, sizeof(snew->sname)-1);
938 snew->sname[sizeof(snew->sname)-1] = '\0';
939 snew->styp = st;
940 snew->priv = NULL;
941 snew->nfargs = n;
942 for (n = 0; n < snew->nfargs; n++)
943 if (fscanf(fp, "%lf", &snew->farg[n]) != 1) {
944 fprintf(stderr, "%s: error reading arguments for '%s'\n",
945 progname, oname);
946 exit(1);
947 }
948 switch (st) {
949 case ST_RING:
950 if (snew->nfargs != 8)
951 goto badcount;
952 VCOPY(snew->snrm, snew->farg+3);
953 if (normalize(snew->snrm) == 0)
954 goto badnorm;
955 if (snew->farg[7] < snew->farg[6]) {
956 double t = snew->farg[7];
957 snew->farg[7] = snew->farg[6];
958 snew->farg[6] = t;
959 }
960 snew->area = PI*(snew->farg[7]*snew->farg[7] -
961 snew->farg[6]*snew->farg[6]);
962 break;
963 case ST_POLY:
964 if (snew->nfargs < 9 || snew->nfargs % 3)
965 goto badcount;
966 finish_polygon(snew);
967 break;
968 case ST_SOURCE:
969 if (snew->nfargs != 4)
970 goto badcount;
971 VCOPY(snew->snrm, snew->farg);
972 if (normalize(snew->snrm) == 0)
973 goto badnorm;
974 snew->area = 2.*PI*(1. - cos((PI/180./2.)*snew->farg[3]));
975 break;
976 }
977 if (snew->area <= FTINY) {
978 fprintf(stderr, "%s: warning - zero area for surface '%s'\n",
979 progname, oname);
980 free(snew);
981 return;
982 }
983 VSUM(curparams.nrm, curparams.nrm, snew->snrm, snew->area);
984 snew->next = curparams.slist;
985 curparams.slist = snew;
986 curparams.nsurfs++;
987 return;
988 badcount:
989 fprintf(stderr, "%s: bad argument count for surface '%s'\n",
990 progname, oname);
991 exit(1);
992 badnorm:
993 fprintf(stderr, "%s: bad orientation for surface '%s'\n",
994 progname, oname);
995 exit(1);
996 }
997
998 /* Parse a receiver object (look for modifiers to add to rcontrib command) */
999 static int
1000 add_recv_object(FILE *fp)
1001 {
1002 int st;
1003 char thismod[128], otype[32], oname[128];
1004 int n;
1005
1006 if (fscanf(fp, "%s %s %s", thismod, otype, oname) != 3)
1007 return(0); /* must have hit EOF! */
1008 if (!strcmp(otype, "alias")) {
1009 fscanf(fp, "%*s"); /* skip alias */
1010 return(0);
1011 }
1012 /* is it a new receiver? */
1013 if ((st = surf_type(otype)) != ST_NONE) {
1014 if (curparams.slist != NULL && (st == ST_SOURCE) ^
1015 (curparams.slist->styp == ST_SOURCE)) {
1016 fputs(progname, stderr);
1017 fputs(": cannot mix source/non-source receivers!\n", stderr);
1018 return(-1);
1019 }
1020 if (strcmp(thismod, curmod)) {
1021 if (curmod[0]) { /* output last receiver? */
1022 finish_receiver();
1023 clear_params(&curparams, 1);
1024 }
1025 strcpy(curmod, thismod);
1026 }
1027 add_surface(st, oname, fp); /* read & store surface */
1028 return(1);
1029 }
1030 /* else skip arguments */
1031 if (!fscanf(fp, "%d", &n)) return;
1032 while (n-- > 0) fscanf(fp, "%*s");
1033 if (!fscanf(fp, "%d", &n)) return;
1034 while (n-- > 0) fscanf(fp, "%*d");
1035 if (!fscanf(fp, "%d", &n)) return;
1036 while (n-- > 0) fscanf(fp, "%*f");
1037 return(0);
1038 }
1039
1040 /* Parse a sender object */
1041 static int
1042 add_send_object(FILE *fp)
1043 {
1044 int st;
1045 char otype[32], oname[128];
1046 int n;
1047
1048 if (fscanf(fp, "%*s %s %s", otype, oname) != 2)
1049 return(0); /* must have hit EOF! */
1050 if (!strcmp(otype, "alias")) {
1051 fscanf(fp, "%*s"); /* skip alias */
1052 return(0);
1053 }
1054 /* is it a new surface? */
1055 if ((st = surf_type(otype)) != ST_NONE) {
1056 if (st == ST_SOURCE) {
1057 fputs(progname, stderr);
1058 fputs(": cannot use source as a sender!\n", stderr);
1059 return(-1);
1060 }
1061 add_surface(st, oname, fp); /* read & store surface */
1062 return(0);
1063 }
1064 /* else skip arguments */
1065 if (!fscanf(fp, "%d", &n)) return;
1066 while (n-- > 0) fscanf(fp, "%*s");
1067 if (!fscanf(fp, "%d", &n)) return;
1068 while (n-- > 0) fscanf(fp, "%*d");
1069 if (!fscanf(fp, "%d", &n)) return;
1070 while (n-- > 0) fscanf(fp, "%*f");
1071 return(0);
1072 }
1073
1074 /* Load a Radiance scene using the given callback function for objects */
1075 static int
1076 load_scene(const char *inspec, int (*ocb)(FILE *))
1077 {
1078 int rv = 0;
1079 char inpbuf[1024];
1080 FILE *fp;
1081 int c;
1082
1083 if (*inspec == '!')
1084 fp = popen(inspec+1, "r");
1085 else
1086 fp = fopen(inspec, "r");
1087 if (fp == NULL) {
1088 fprintf(stderr, "%s: cannot load '%s'\n", progname, inspec);
1089 return(-1);
1090 }
1091 while ((c = getc(fp)) != EOF) { /* load receiver data */
1092 if (isspace(c)) /* skip leading white space */
1093 continue;
1094 if (c == '!') { /* read from a new command */
1095 inpbuf[0] = c;
1096 if (fgetline(inpbuf+1, sizeof(inpbuf)-1, fp) != NULL) {
1097 if ((c = load_scene(inpbuf, ocb)) < 0)
1098 return(c);
1099 rv += c;
1100 }
1101 continue;
1102 }
1103 if (c == '#') { /* parameters/comment */
1104 if (fscanf(fp, "%s", inpbuf) == 1 &&
1105 !strcmp(inpbuf, PARAMSTART)) {
1106 if (fgets(inpbuf, sizeof(inpbuf), fp) != NULL)
1107 parse_params(inpbuf);
1108 continue;
1109 }
1110 while ((c = getc(fp)) != EOF && c != '\n');
1111 ; /* else skipping comment */
1112 continue;
1113 }
1114 ungetc(c, fp); /* else check object for receiver */
1115 c = (*ocb)(fp);
1116 if (c < 0)
1117 return(c);
1118 rv += c;
1119 }
1120 /* close our input stream */
1121 c = (*inspec == '!') ? pclose(fp) : fclose(fp);
1122 if (c != 0) {
1123 fprintf(stderr, "%s: error loading '%s'\n", progname, inspec);
1124 return(-1);
1125 }
1126 return(rv);
1127 }
1128
1129 /* Get command arguments and run program */
1130 int
1131 main(int argc, char *argv[])
1132 {
1133 char fmtopt[6] = "-faa"; /* default output is ASCII */
1134 char *sendfn;
1135 char sampcntbuf[32], nsbinbuf[32];
1136 FILE *rcfp;
1137 int nsbins;
1138 int a, i;
1139 /* screen rcontrib options */
1140 progname = argv[0];
1141 for (a = 1; a < argc-2 && argv[a][0] == '-'; a++) {
1142 int na = 1; /* !! Keep consistent !! */
1143 switch (argv[a][1]) {
1144 case 'v': /* verbose mode */
1145 verbose = !verbose;
1146 na = 0;
1147 continue;
1148 case 'f': /* special case for -fo, -ff, etc. */
1149 switch (argv[a][2]) {
1150 case '\0': /* cal file */
1151 goto userr;
1152 case 'o': /* force output */
1153 goto userr;
1154 case 'a': /* output format */
1155 case 'f':
1156 case 'd':
1157 case 'c':
1158 if (!(fmtopt[4] = argv[a][3]))
1159 fmtopt[4] = argv[a][2];
1160 fmtopt[3] = argv[a][2];
1161 na = 0;
1162 continue; /* will pass later */
1163 default:
1164 goto userr;
1165 }
1166 break;
1167 case 'c': /* number of samples */
1168 sampcnt = atoi(argv[a+1]);
1169 if (sampcnt <= 0)
1170 goto userr;
1171 na = 0; /* we re-add this later */
1172 continue;
1173 case 'V': /* options without arguments */
1174 case 'w':
1175 case 'u':
1176 case 'i':
1177 case 'h':
1178 break;
1179 case 'n': /* options with 1 argument */
1180 case 's':
1181 case 'o':
1182 na = 2;
1183 break;
1184 case 'b': /* special case */
1185 if (argv[a][2] != 'v') goto userr;
1186 break;
1187 case 'l': /* special case */
1188 if (argv[a][2] == 'd') goto userr;
1189 na = 2;
1190 break;
1191 case 'd': /* special case */
1192 if (argv[a][2] != 'v') na = 2;
1193 break;
1194 case 'a': /* special case */
1195 na = (argv[a][2] == 'v') ? 4 : 2;
1196 break;
1197 case 'm': /* special case */
1198 if (!argv[a][2]) goto userr;
1199 na = (argv[a][2] == 'e') | (argv[a][2] == 'a') ? 4 : 2;
1200 break;
1201 case '\0': /* pass-through mode */
1202 goto done_opts;
1203 default: /* anything else is verbotten */
1204 goto userr;
1205 }
1206 if (na <= 0) continue;
1207 CHECKARGC(na); /* pass on option */
1208 rcarg[nrcargs++] = argv[a];
1209 while (--na) /* + arguments if any */
1210 rcarg[nrcargs++] = argv[++a];
1211 }
1212 done_opts:
1213 if (a > argc-2)
1214 goto userr; /* check at end of options */
1215 sendfn = argv[a++]; /* assign sender & receiver inputs */
1216 if (sendfn[0] == '-') { /* user wants pass-through mode? */
1217 if (sendfn[1]) goto userr;
1218 sendfn = NULL;
1219 if (sampcnt <= 0) sampcnt = 1;
1220 } else { /* else FVECT determines input format */
1221 fmtopt[3] = (sizeof(RREAL)==sizeof(double)) ? 'd' : 'f';
1222 if (sampcnt <= 0) sampcnt = 10000;
1223 }
1224 sprintf(sampcntbuf, "%d", sampcnt);
1225 CHECKARGC(3); /* add our format & sample count */
1226 rcarg[nrcargs++] = fmtopt;
1227 rcarg[nrcargs++] = "-c";
1228 rcarg[nrcargs++] = sampcntbuf;
1229 /* add receiver arguments to rcontrib */
1230 if (load_scene(argv[a], add_recv_object) < 0)
1231 return(1);
1232 finish_receiver();
1233 if (sendfn == NULL) { /* pass-through mode? */
1234 CHECKARGC(1); /* add octree */
1235 rcarg[nrcargs++] = oconv_command(argc-a, argv+a);
1236 rcarg[nrcargs] = NULL;
1237 return(my_exec(rcarg)); /* rcontrib does everything */
1238 }
1239 clear_params(&curparams, 0); /* else load sender surface & params */
1240 if (load_scene(sendfn, add_send_object) < 0)
1241 return(1);
1242 if ((nsbins = prepare_sampler()) <= 0)
1243 return(1);
1244 CHECKARGC(3); /* add row count and octree */
1245 rcarg[nrcargs++] = "-y";
1246 sprintf(nsbinbuf, "%d", nsbins);
1247 rcarg[nrcargs++] = nsbinbuf;
1248 rcarg[nrcargs++] = oconv_command(argc-a, argv+a);
1249 rcarg[nrcargs] = NULL;
1250 /* open pipe to rcontrib process */
1251 if ((rcfp = popen_arglist(rcarg, "w")) == NULL)
1252 return(1);
1253 SET_FILE_BINARY(rcfp);
1254 #ifdef getc_unlocked
1255 flockfile(rcfp);
1256 #endif
1257 if (verbose) {
1258 fprintf(stderr, "%s: sampling %d directions", progname, nsbins);
1259 if (curparams.nsurfs > 1)
1260 fprintf(stderr, " (%d surfaces)\n", curparams.nsurfs);
1261 else
1262 fputc('\n', stderr);
1263 }
1264 for (i = 0; i < nsbins; i++) /* send rcontrib ray samples */
1265 if (!(*curparams.sample_basis)(&curparams, i, rcfp))
1266 return(1);
1267 return(pclose(rcfp) == 0); /* all finished! */
1268 userr:
1269 if (a < argc-2)
1270 fprintf(stderr, "%s: unsupported option '%s'", progname, argv[a]);
1271 fprintf(stderr, "Usage: %s [rcontrib options] sender.rad receiver.rad [system.rad ..]\n",
1272 progname);
1273 return(1);
1274 }