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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.2 static const char RCSid[] = "$Id: rfluxmtx.c,v 2.1 2014/07/19 01:19:49 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     char *binv = NULL;
377     char *binf = NULL;
378     char *nbins = NULL;
379     char sbuf[256];
380 greg 2.1
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 greg 2.2 binv = "0";
409 greg 2.1 } else if (tolower(curparams.hemis[0]) == 's' &&
410     tolower(curparams.hemis[1]) == 'c') {
411 greg 2.2 /* assign parameters */
412 greg 2.1 if (curparams.hsiz <= 1) {
413     fputs(progname, stderr);
414     fputs(": missing size for Shirley-Chiu sampling!\n", stderr);
415     exit(1);
416     }
417 greg 2.2 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 greg 2.1 } else if ((tolower(curparams.hemis[0]) == 'r') |
428     (tolower(curparams.hemis[0]) == 't')) {
429 greg 2.2 calfn = reinhfn; reinhfn = NULL;
430     /* XXX Need to set number of divisions */
431     binf = "rhbin";
432     nbins = "Nrhbins";
433 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
434     !curparams.hemis[1] |
435     (tolower(curparams.hemis[1]) == 'f') |
436     (curparams.hemis[1] == '1')) {
437 greg 2.2 calfn = kfullfn; kfullfn = NULL;
438     binf = "kbin";
439     nbins = "Nkbins";
440 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
441     (tolower(curparams.hemis[1]) == 'h') |
442     (curparams.hemis[1] == '2')) {
443 greg 2.2 calfn = khalffn; khalffn = NULL;
444     binf = "khbin";
445     nbins = "Nkhbins";
446 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
447     (tolower(curparams.hemis[1]) == 'q') |
448     (curparams.hemis[1] == '4')) {
449 greg 2.2 calfn = kquarterfn; kquarterfn = NULL;
450     binf = "kqbin";
451     nbins = "Nkqbins";
452 greg 2.1 } else {
453     fprintf(stderr, "%s: unrecognized hemisphere sampling: h=%s\n",
454     progname, curparams.hemis);
455     exit(1);
456     }
457 greg 2.2 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 greg 2.1 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     }