ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.3
Committed: Tue Jul 22 02:12:48 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +19 -11 lines
Log Message:
Added Reinhart sky cal file

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.3 static const char RCSid[] = "$Id: rfluxmtx.c,v 2.2 2014/07/21 15:59:47 greg Exp $";
3 greg 2.1 #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 greg 2.2 #define SPECIALS " \t\"$*?"
27 greg 2.1 #define QUOTCHAR '"'
28     #else
29 greg 2.2 #define SPECIALS " \t\n'\"()${}*?[];|&"
30 greg 2.1 #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 greg 2.2 int sampcnt = 0; /* sample count (0==unset) */
49 greg 2.1
50 greg 2.2 char *reinhfn = "reinhartb.cal";
51     char *shirchiufn = "disk2square.cal";
52 greg 2.1 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 greg 2.2 short nfargs; /* number of real arguments */
72     double farg[1]; /* real values (extends struct) */
73 greg 2.1 } 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 greg 2.2 char *calfn = NULL;
376 greg 2.3 char *params = NULL;
377 greg 2.2 char *binv = NULL;
378     char *binf = NULL;
379     char *nbins = NULL;
380     char sbuf[256];
381 greg 2.1
382     if (!curmod[0]) {
383     fputs(progname, stderr);
384     fputs(": missing receiver surface!\n", stderr);
385     exit(1);
386     }
387     if (curparams.outfn != NULL) { /* add output file spec. */
388     CHECKARGC(2);
389     rcarg[nrcargs++] = "-o";
390     rcarg[nrcargs++] = curparams.outfn;
391     }
392     /* add bin specification */
393     if (!curparams.hemis[0]) {
394     fputs(progname, stderr);
395     fputs(": missing hemisphere sampling type!\n", stderr);
396     exit(1);
397     }
398     if (normalize(curparams.nrm) == 0) {
399     fputs(progname, stderr);
400     fputs(": undefined normal for hemisphere sampling\n", stderr);
401     exit(1);
402     }
403     if (normalize(curparams.vup) == 0)
404     if (fabs(curparams.nrm[2]) < .7)
405     curparams.vup[2] = 1;
406     else
407     curparams.vup[1] = 1;
408     if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1') {
409 greg 2.3 binv = "0"; /* uniform sampling -- one bin */
410 greg 2.1 } else if (tolower(curparams.hemis[0]) == 's' &&
411     tolower(curparams.hemis[1]) == 'c') {
412 greg 2.2 /* assign parameters */
413 greg 2.1 if (curparams.hsiz <= 1) {
414     fputs(progname, stderr);
415     fputs(": missing size for Shirley-Chiu sampling!\n", stderr);
416     exit(1);
417     }
418 greg 2.3 calfn = shirchiufn; shirchiufn = NULL;
419 greg 2.2 sprintf(sbuf, "SCdim=%d,Nx=%g,Ny=%g,Nz=%g,Ux=%g,Uy=%g,Uz=%g",
420     curparams.hsiz,
421     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
422     curparams.vup[0], curparams.vup[1], curparams.vup[2]);
423 greg 2.3 params = savqstr(sbuf);
424 greg 2.2 binv = "scbin";
425     nbins = "SCdim*SCdim";
426 greg 2.1 } else if ((tolower(curparams.hemis[0]) == 'r') |
427     (tolower(curparams.hemis[0]) == 't')) {
428 greg 2.2 calfn = reinhfn; reinhfn = NULL;
429 greg 2.3 sprintf(sbuf, "MF=%d,Nx=%g,Ny=%g,Nz=%g,Ux=%g,Uy=%g,Uz=%g",
430     curparams.hsiz,
431     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
432     curparams.vup[0], curparams.vup[1], curparams.vup[2]);
433     params = savqstr(sbuf);
434     binv = "rbin";
435     nbins = "Nrbins";
436 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
437     !curparams.hemis[1] |
438     (tolower(curparams.hemis[1]) == 'f') |
439     (curparams.hemis[1] == '1')) {
440 greg 2.2 calfn = kfullfn; kfullfn = NULL;
441     binf = "kbin";
442     nbins = "Nkbins";
443 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
444     (tolower(curparams.hemis[1]) == 'h') |
445     (curparams.hemis[1] == '2')) {
446 greg 2.2 calfn = khalffn; khalffn = NULL;
447     binf = "khbin";
448     nbins = "Nkhbins";
449 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
450     (tolower(curparams.hemis[1]) == 'q') |
451     (curparams.hemis[1] == '4')) {
452 greg 2.2 calfn = kquarterfn; kquarterfn = NULL;
453     binf = "kqbin";
454     nbins = "Nkqbins";
455 greg 2.1 } else {
456     fprintf(stderr, "%s: unrecognized hemisphere sampling: h=%s\n",
457     progname, curparams.hemis);
458     exit(1);
459     }
460 greg 2.2 if (calfn != NULL) { /* add cal file if needed */
461     CHECKARGC(2);
462     rcarg[nrcargs++] = "-f";
463     rcarg[nrcargs++] = calfn;
464     }
465 greg 2.3 if (params != NULL) { /* parameters _after_ cal file */
466     CHECKARGC(2);
467     rcarg[nrcargs++] = "-p";
468     rcarg[nrcargs++] = params;
469     }
470 greg 2.2 if (nbins != NULL) { /* add #bins if set */
471     CHECKARGC(2);
472     rcarg[nrcargs++] = "-bn";
473     rcarg[nrcargs++] = nbins;
474     }
475 greg 2.3 if (binv != NULL) {
476 greg 2.2 CHECKARGC(2); /* assign bin variable */
477     rcarg[nrcargs++] = "-b";
478     rcarg[nrcargs++] = binv;
479     } else if (binf != NULL) {
480     CHECKARGC(2); /* assign bin function */
481 greg 2.3 rcarg[nrcargs++] = "-b";
482 greg 2.2 sprintf(sbuf, "%s(%g,%g,%g,%g,%g,%g)", binf,
483     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
484     curparams.vup[0], curparams.vup[1], curparams.vup[2]);
485     rcarg[nrcargs++] = savqstr(sbuf);
486     }
487     CHECKARGC(2); /* modifier argument goes last */
488 greg 2.1 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     }