ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/rt/rxfluxmtx.cpp
Revision: 2.1
Committed: Wed Oct 22 02:38:27 2025 UTC (2 weeks, 1 day ago) by greg
Branch: MAIN
Log Message:
feat(rxfluxmtx): Created new tool with more x's in its name than any other

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2     static const char RCSid[] = "$Id: rfluxmtx.c,v 2.61 2025/10/17 17:24:47 greg Exp $";
3     #endif
4     /*
5     * Calculate flux transfer matrix or matrices using RcontribSimulManager class
6     *
7     * Much of this is cribbed from util/rfluxmtx.c
8     * Implementation here improves efficiency and overcomes receiver limits
9     */
10    
11     #include "copyright.h"
12    
13     #include <ctype.h>
14     #include <signal.h>
15     #include "RcontribSimulManager.h"
16     #include "bsdf.h"
17     #include "bsdf_m.h"
18     #include "func.h"
19     #include "random.h"
20     #include "triangulate.h"
21     #include "view.h"
22    
23     int verby = 0; /* verbose mode? */
24    
25     const char *reinhfn = "reinhartb.cal";
26     const char *shirchiufn = "disk2square.cal";
27     const char *kfullfn = "klems_full.cal";
28     const char *khalffn = "klems_half.cal";
29     const char *kquarterfn = "klems_quarter.cal";
30     const char *ciefn = "cieskyscan.cal";
31    
32     const char *sigerr[NSIG]; /* signal error messages */
33    
34     int inpfmt = 'a'; /* input format */
35     int outfmt = 'f'; /* output format */
36    
37     RcontribSimulManager myRCmanager; // global rcontrib simulation manager
38    
39     #define PARAMSNAME "rfluxmtx" /* string indicating parameters */
40    
41     /* surface type IDs */
42     #define ST_NONE 0
43     #define ST_POLY 1
44     #define ST_RING 2
45     #define ST_SOURCE 3
46     /* surface structure */
47     struct SURF {
48     SURF *next; /* next surface in list */
49     void *priv; /* private data (malloc'ed) */
50     char sname[32]; /* surface name */
51     FVECT snrm; /* surface normal */
52     double area; /* surface area / proj. solid angle */
53     short styp; /* surface type */
54     short nfargs; /* number of real arguments */
55     double farg[1]; /* real values (extends struct) */
56     };
57     /* triangle */
58     struct PTRI {
59     double afrac; /* fraction of total area */
60     short vndx[3]; /* vertex indices */
61     };
62     /* triangulated polygon */
63     struct POLYTRIS {
64     FVECT uva[2]; /* tangent axes */
65     int ntris; /* number of triangles */
66     PTRI tri[1]; /* triangle array (extends struct) */
67     };
68     /* sender/receiver parameters */
69     struct PARAMS {
70     char sign; /* '-' for axis reversal */
71     char hemis[31]; /* hemispherical sampling spec. */
72     int hsiz; /* hemisphere basis size */
73     int nsurfs; /* number of surfaces */
74     SURF *slist; /* list of surfaces */
75     FVECT vup; /* up vector (zero if unset) */
76     FVECT nrm; /* average normal direction */
77     FVECT udir, vdir; /* tangent axes */
78     char outfn[256]; /* output file name (receiver) */
79     int (*sample_basis)(PARAMS *p, int b, FVECT orig_dir[]);
80     };
81    
82     PARAMS curparams;
83     char curmod[MAXSTR];
84     char newparams[1024];
85    
86     typedef int SURFSAMP(FVECT, SURF *, double);
87    
88     SURFSAMP ssamp_bad, ssamp_poly, ssamp_ring;
89    
90     SURFSAMP *orig_in_surf[4] = {
91     ssamp_bad, ssamp_poly, ssamp_ring, ssamp_bad
92     };
93    
94     VIEW ourview = STDVIEW; // view parameters (if set)
95     double dstrpix = 0; // pixel jitter
96     double dblur = 0; // depth-of-field
97    
98     /* Clear parameter set */
99     void
100     clear_params(PARAMS *p, bool reset_only = true)
101     {
102     while (p->slist != NULL) {
103     SURF *sdel = p->slist;
104     p->slist = sdel->next;
105     if (sdel->priv != NULL)
106     free(sdel->priv);
107     free(sdel);
108     }
109     if (reset_only) {
110     p->nsurfs = 0;
111     p->slist = NULL;
112     memset(p->nrm, 0, sizeof(FVECT));
113     memset(p->vup, 0, sizeof(FVECT));
114     } else
115     memset(p, 0, sizeof(PARAMS));
116     }
117    
118     /* Get surface type from name */
119     int
120     surf_type(const char *otype)
121     {
122     if (!strcmp(otype, "polygon"))
123     return(ST_POLY);
124     if (!strcmp(otype, "ring"))
125     return(ST_RING);
126     if (!strcmp(otype, "source"))
127     return(ST_SOURCE);
128     return(ST_NONE);
129     }
130    
131     /* Add arguments to oconv command */
132     char *
133     oconv_command(int ac, char *av[])
134     {
135     static char oconvbuf[4096] = "!oconv -f ";
136     char *cp = oconvbuf + 10;
137     char *recv = *av++;
138    
139     if (ac-- <= 0)
140     return(NULL);
141     if (erract[WARNING].pf == NULL) { /* warnings off? */
142     strcpy(cp, "-w ");
143     cp += 3;
144     }
145     while (ac-- > 0) { /* copy each argument */
146     int len = strlen(*av);
147     if (cp+len+4 >= oconvbuf+sizeof(oconvbuf))
148     goto overrun;
149     if (matchany(*av, SPECIALS)) {
150     *cp++ = QUOTCHAR;
151     strcpy(cp, *av++);
152     cp += len;
153     *cp++ = QUOTCHAR;
154     } else {
155     strcpy(cp, *av++);
156     cp += len;
157     }
158     *cp++ = ' ';
159     }
160     /* receiver goes last */
161     if (matchany(recv, SPECIALS)) {
162     *cp++ = QUOTCHAR;
163     while (*recv) {
164     if (cp >= oconvbuf+(sizeof(oconvbuf)-3))
165     goto overrun;
166     *cp++ = *recv++;
167     }
168     *cp++ = QUOTCHAR;
169     *cp = '\0';
170     } else
171     strcpy(cp, recv);
172     return(oconvbuf);
173     overrun:
174     error(USER, "too many file arguments!");
175     return(NULL); /* pro forma return */
176     }
177    
178     /* Get normalized direction vector from string specification */
179     int
180     get_direction(FVECT dv, const char *s)
181     {
182     int sign = 1;
183     int axis = 0;
184    
185     memset(dv, 0, sizeof(FVECT));
186     nextchar:
187     switch (*s) {
188     case '+':
189     ++s;
190     goto nextchar;
191     case '-':
192     sign = -sign;
193     ++s;
194     goto nextchar;
195     case 'z':
196     case 'Z':
197     ++axis;
198     case 'y':
199     case 'Y':
200     ++axis;
201     case 'x':
202     case 'X':
203     dv[axis] = sign;
204     return(!s[1] | isspace(s[1]));
205     default:
206     break;
207     }
208     #ifdef SMLFLT
209     if (sscanf(s, "%f,%f,%f", &dv[0], &dv[1], &dv[2]) != 3)
210     #else
211     if (sscanf(s, "%lf,%lf,%lf", &dv[0], &dv[1], &dv[2]) != 3)
212     #endif
213     return(0);
214     dv[0] *= (RREAL)sign;
215     return(normalize(dv) > 0);
216     }
217    
218     /* Parse program parameters (directives) */
219     int
220     parse_params(PARAMS *p, char *pargs)
221     {
222     char *cp = pargs;
223     int nparams = 0;
224     int quot;
225     int i;
226    
227     for ( ; ; ) {
228     switch (*cp++) {
229     case 'h':
230     if (*cp++ != '=')
231     break;
232     if ((*cp == '+') | (*cp == '-'))
233     p->sign = *cp++;
234     else
235     p->sign = '+';
236     p->hsiz = 0;
237     i = 0;
238     while (*cp && !isspace(*cp)) {
239     if (isdigit(*cp))
240     p->hsiz = 10*p->hsiz + *cp - '0';
241     p->hemis[i++] = *cp++;
242     }
243     if (!i)
244     break;
245     p->hemis[i] = '\0';
246     p->hsiz += !p->hsiz;
247     ++nparams;
248     continue;
249     case 'u':
250     if (*cp++ != '=')
251     break;
252     if (!get_direction(p->vup, cp))
253     break;
254     while (*cp && !isspace(*cp++))
255     ;
256     ++nparams;
257     continue;
258     case 'o':
259     if (*cp++ != '=')
260     break;
261     quot = 0;
262     if ((*cp == '"') | (*cp == '\''))
263     quot = *cp++;
264     i = 0;
265     while (*cp && (quot ? (*cp != quot) : !isspace(*cp))) {
266     i++; cp++;
267     }
268     if (!i)
269     break;
270     if (!*cp) {
271     if (quot)
272     break;
273     cp[1] = '\0';
274     }
275     *cp = '\0';
276     strlcpy(p->outfn, cp-i, sizeof(p->outfn));
277     *cp++ = quot ? quot : ' ';
278     ++nparams;
279     continue;
280     case ' ':
281     case '\t':
282     case '\r':
283     case '\n':
284     continue;
285     case '\0':
286     return(nparams);
287     default:
288     break;
289     }
290     break;
291     }
292     sprintf(errmsg, "bad parameter string: %s", pargs);
293     error(USER, errmsg);
294     return(-1); /* pro forma return */
295     }
296    
297     /* Add receiver modifier and associated parameters */
298     void
299     finish_receiver()
300     {
301     bool uniform = false;
302     const char *calfn = NULL;
303     char binv[64] = "";
304     char params[128] = "";
305     const char *binf = NULL;
306     const char *nbins = NULL;
307     /* check arguments */
308     if (!curmod[0])
309     error(USER, "missing receiver surface");
310     if (!curparams.outfn[0])
311     error(USER, "missing output file for receiver");
312     if (!curparams.hemis[0])
313     error(USER, "missing hemisphere sampling type");
314     if (normalize(curparams.nrm) == 0)
315     error(USER, "undefined normal for hemisphere sampling");
316     if (normalize(curparams.vup) == 0) {
317     if (fabs(curparams.nrm[2]) < .7)
318     curparams.vup[2] = 1;
319     else
320     curparams.vup[1] = 1;
321     }
322     /* determine sample type/bin */
323     if ((tolower(curparams.hemis[0]) == 'u') | (curparams.hemis[0] == '1')) {
324     sprintf(binv, "if(-Dx*%g-Dy*%g-Dz*%g,0,-1)",
325     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2]);
326     uniform = true; /* uniform sampling -- one bin */
327     } else if (tolower(curparams.hemis[0]) == 's' &&
328     tolower(curparams.hemis[1]) == 'c') {
329     /* assign parameters */
330     if (curparams.hsiz <= 1)
331     error(USER, "missing size for Shirley-Chiu sampling");
332     calfn = shirchiufn; shirchiufn = NULL;
333     sprintf(params, "SCdim=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g,RHS=%c1",
334     curparams.hsiz,
335     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
336     curparams.vup[0], curparams.vup[1], curparams.vup[2],
337     curparams.sign);
338     strcpy(binv, "scbin");
339     nbins = "SCdim*SCdim";
340     } else if ((tolower(curparams.hemis[0]) == 'r') |
341     (tolower(curparams.hemis[0]) == 't')) {
342     calfn = reinhfn; reinhfn = NULL;
343     sprintf(params, "MF=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g,RHS=%c1",
344     curparams.hsiz,
345     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
346     curparams.vup[0], curparams.vup[1], curparams.vup[2],
347     curparams.sign);
348     strcpy(binv, "rbin");
349     nbins = "Nrbins";
350     } else if (tolower(curparams.hemis[0]) == 'k' &&
351     !curparams.hemis[1] |
352     (tolower(curparams.hemis[1]) == 'f') |
353     (curparams.hemis[1] == '1')) {
354     calfn = kfullfn; kfullfn = NULL;
355     binf = "kbin";
356     nbins = "Nkbins";
357     } else if (tolower(curparams.hemis[0]) == 'k' &&
358     (tolower(curparams.hemis[1]) == 'h') |
359     (curparams.hemis[1] == '2')) {
360     calfn = khalffn; khalffn = NULL;
361     binf = "khbin";
362     nbins = "Nkhbins";
363     } else if (tolower(curparams.hemis[0]) == 'k' &&
364     (tolower(curparams.hemis[1]) == 'q') |
365     (curparams.hemis[1] == '4')) {
366     calfn = kquarterfn; kquarterfn = NULL;
367     binf = "kqbin";
368     nbins = "Nkqbins";
369     } else if (!strcasecmp(curparams.hemis, "cie")) {
370     calfn = ciefn; ciefn = NULL;
371     sprintf(params, "rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g,RHS=%c1",
372     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
373     curparams.vup[0], curparams.vup[1], curparams.vup[2],
374     curparams.sign);
375     strcpy(binv, "cbin");
376     nbins = "Ncbins";
377     } else {
378     sprintf(errmsg, "unrecognized hemisphere sampling: h=%s",
379     curparams.hemis);
380     error(USER, errmsg);
381     }
382     if (tolower(curparams.hemis[0]) == 'k') {
383     sprintf(params, "RHS=%c1", curparams.sign);
384     }
385     if (!uniform)
386     for (SURF *sp = curparams.slist; sp != NULL; sp = sp->next)
387     if (sp->styp == ST_SOURCE && fabs(sp->area - PI) > 1e-3) {
388     sprintf(errmsg, "source '%s' must be 180-degrees",
389     sp->sname);
390     error(USER, errmsg);
391     }
392     if (binf != NULL) { /* complete bin function? */
393     sprintf(binv, "%s(%g,%g,%g,%g,%g,%g)", binf,
394     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
395     curparams.vup[0], curparams.vup[1], curparams.vup[2]);
396     }
397     if (calfn != NULL) /* load cal file if needed */
398     loadfunc(const_cast<char *>(calfn));
399     if (params[0]) /* set parameters if any */
400     set_eparams(params);
401     // add modifier to rcontrib object
402     if (!myRCmanager.AddModifier(curmod, curparams.outfn, params, binv,
403     nbins==NULL ? 1 : int(eval(const_cast<char *>(nbins)) + .5)))
404     error(INTERNAL, "AddModifier() call failed");
405    
406     clear_params(&curparams, true); /* clear for next receiver */
407     curmod[0] = '\0';
408     }
409    
410     /* Make randomly oriented tangent plane axes for given normal direction */
411     void
412     make_axes(FVECT uva[2], const FVECT nrm)
413     {
414     int i;
415    
416     if (!getperpendicular(uva[0], nrm, 1))
417     error(USER, "bad surface normal in make_axes()");
418     VCROSS(uva[1], nrm, uva[0]);
419     }
420    
421     /* Illegal sender surfaces end up here */
422     int
423     ssamp_bad(FVECT orig, SURF *sp, double x)
424     {
425     sprintf(errmsg, "illegal sender surface '%s'", sp->sname);
426     error(INTERNAL, errmsg);
427     return(0);
428     }
429    
430     /* Generate origin on ring surface from uniform random variable */
431     int
432     ssamp_ring(FVECT orig, SURF *sp, double x)
433     {
434     FVECT *uva = (FVECT *)sp->priv;
435     double samp2[2];
436     double uv[2];
437     int i;
438    
439     if (uva == NULL) { /* need tangent axes */
440     uva = (FVECT *)malloc(sizeof(FVECT)*2);
441     if (uva == NULL) {
442     error(SYSTEM, "out of memory in ssamp_ring()");
443     return(0);
444     }
445     make_axes(uva, sp->snrm);
446     sp->priv = uva;
447     }
448     SDmultiSamp(samp2, 2, x);
449     samp2[0] = sqrt(samp2[0]*sp->area*(1./PI) + sp->farg[6]*sp->farg[6]);
450     samp2[1] *= 2.*PI;
451     uv[0] = samp2[0]*tcos(samp2[1]);
452     uv[1] = samp2[0]*tsin(samp2[1]);
453     for (i = 3; i--; )
454     orig[i] = sp->farg[i] + uv[0]*uva[0][i] + uv[1]*uva[1][i];
455     return(1);
456     }
457    
458     /* Add triangle to polygon's list (call-back function) */
459     int
460     add_triangle(const Vert2_list *tp, int a, int b, int c)
461     {
462     POLYTRIS *ptp = (POLYTRIS *)tp->p;
463     PTRI *trip = ptp->tri + ptp->ntris++;
464    
465     trip->vndx[0] = a;
466     trip->vndx[1] = b;
467     trip->vndx[2] = c;
468     return(1);
469     }
470    
471     /* Generate origin on polygon surface from uniform random variable */
472     int
473     ssamp_poly(FVECT orig, SURF *sp, double x)
474     {
475     POLYTRIS *ptp = (POLYTRIS *)sp->priv;
476     double samp2[2];
477     double *v0, *v1, *v2;
478     int i;
479    
480     if (ptp == NULL) { /* need to triangulate */
481     ptp = (POLYTRIS *)malloc(sizeof(POLYTRIS) +
482     sizeof(PTRI)*(sp->nfargs/3 - 3));
483     if (ptp == NULL)
484     goto memerr;
485     if (sp->nfargs == 3) { /* simple case */
486     ptp->ntris = 1;
487     ptp->tri[0].vndx[0] = 0;
488     ptp->tri[0].vndx[1] = 1;
489     ptp->tri[0].vndx[2] = 2;
490     ptp->tri[0].afrac = 1;
491     } else {
492     Vert2_list *v2l = polyAlloc(sp->nfargs/3);
493     if (v2l == NULL)
494     goto memerr;
495     make_axes(ptp->uva, sp->snrm);
496     for (i = v2l->nv; i--; ) {
497     v2l->v[i].mX = DOT(sp->farg+3*i, ptp->uva[0]);
498     v2l->v[i].mY = DOT(sp->farg+3*i, ptp->uva[1]);
499     }
500     ptp->ntris = 0;
501     v2l->p = ptp;
502     if (!polyTriangulate(v2l, add_triangle)) {
503     sprintf(errmsg, "cannot triangulate polygon '%s'",
504     sp->sname);
505     error(USER, errmsg);
506     return(0);
507     }
508     for (i = ptp->ntris; i--; ) {
509     int a = ptp->tri[i].vndx[0];
510     int b = ptp->tri[i].vndx[1];
511     int c = ptp->tri[i].vndx[2];
512     ptp->tri[i].afrac =
513     (v2l->v[a].mX*v2l->v[b].mY -
514     v2l->v[b].mX*v2l->v[a].mY +
515     v2l->v[b].mX*v2l->v[c].mY -
516     v2l->v[c].mX*v2l->v[b].mY +
517     v2l->v[c].mX*v2l->v[a].mY -
518     v2l->v[a].mX*v2l->v[c].mY) /
519     (2.*sp->area);
520     }
521     polyFree(v2l);
522     }
523     sp->priv = ptp;
524     }
525     /* pick triangle by partial area */
526     for (i = 0; i < ptp->ntris-1 && x > ptp->tri[i].afrac; i++)
527     x -= ptp->tri[i].afrac;
528     SDmultiSamp(samp2, 2, x/ptp->tri[i].afrac);
529     samp2[0] *= samp2[1] = sqrt(samp2[1]);
530     samp2[1] = 1. - samp2[1];
531     v0 = sp->farg + 3*ptp->tri[i].vndx[0];
532     v1 = sp->farg + 3*ptp->tri[i].vndx[1];
533     v2 = sp->farg + 3*ptp->tri[i].vndx[2];
534     for (i = 3; i--; )
535     orig[i] = v0[i] + samp2[0]*(v1[i] - v0[i])
536     + samp2[1]*(v2[i] - v0[i]) ;
537     return(1);
538     memerr:
539     error(SYSTEM, "out of memory in ssamp_poly");
540     return(0);
541     }
542    
543     /* Compute sample origin based on projected areas of sender subsurfaces */
544     int
545     sample_origin(PARAMS *p, FVECT orig, const FVECT rdir, double x)
546     {
547     static double *projsa;
548     static int nall;
549     double tarea = 0;
550     int i;
551     SURF *sp;
552     /* special case for lone surface */
553     if (p->nsurfs == 1) {
554     sp = p->slist;
555     if (DOT(sp->snrm, rdir) >= FTINY) {
556     sprintf(errmsg, "sample behind sender '%s'", sp->sname);
557     error(INTERNAL, errmsg);
558     return(0);
559     }
560     return((*orig_in_surf[sp->styp])(orig, sp, x));
561     }
562     if (p->nsurfs > nall) { /* (re)allocate surface area cache */
563     if (projsa) free(projsa);
564     projsa = (double *)malloc(sizeof(double)*p->nsurfs);
565     if (projsa == NULL)
566     error(SYSTEM, "out of memory in sample_origin()");
567     nall = p->nsurfs;
568     }
569     /* compute projected areas */
570     for (i = 0, sp = p->slist; sp != NULL; i++, sp = sp->next) {
571     projsa[i] = -DOT(sp->snrm, rdir) * sp->area;
572     tarea += projsa[i] *= (double)(projsa[i] > 0);
573     }
574     if (tarea < FTINY*FTINY) { /* wrong side of sender? */
575     error(INTERNAL, "sample behind all sender elements!");
576     return(0);
577     }
578     tarea *= x; /* get surface from list */
579     for (i = 0, sp = p->slist; tarea > projsa[i]; sp = sp->next)
580     tarea -= projsa[i++];
581     return((*orig_in_surf[sp->styp])(orig, sp, tarea/projsa[i]));
582     }
583    
584     /* Uniform sample generator */
585     int
586     sample_uniform(PARAMS *p, int b, FVECT orig_dir[])
587     {
588     int n = myRCmanager.accum;
589     double samp3[3];
590     FVECT duvw;
591     int i;
592    
593     if (orig_dir == NULL) /* just requesting number of bins? */
594     return(1);
595    
596     while (n--) { /* stratified hemisphere sampling */
597     SDmultiSamp(samp3, 3, (n+frandom())/myRCmanager.accum);
598     square2disk(duvw, samp3[1], samp3[2]);
599     duvw[2] = -sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
600     for (i = 3; i--; )
601     orig_dir[1][i] = duvw[0]*p->udir[i] +
602     duvw[1]*p->vdir[i] +
603     duvw[2]*p->nrm[i] ;
604     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
605     return(0);
606     orig_dir += 2;
607     }
608     return(1);
609     }
610    
611     /* Shirly-Chiu sample generator */
612     int
613     sample_shirchiu(PARAMS *p, int b, FVECT orig_dir[])
614     {
615     int n = myRCmanager.accum;
616     double samp3[3];
617     FVECT duvw;
618     int i;
619    
620     if (orig_dir == NULL) /* just requesting number of bins? */
621     return(p->hsiz*p->hsiz);
622    
623     while (n--) { /* stratified sampling */
624     SDmultiSamp(samp3, 3, (n+frandom())/myRCmanager.accum);
625     square2disk(duvw, (b/p->hsiz + samp3[1])/curparams.hsiz,
626     (b%p->hsiz + samp3[2])/curparams.hsiz);
627     duvw[2] = sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
628     for (i = 3; i--; )
629     orig_dir[1][i] = -duvw[0]*p->udir[i] -
630     duvw[1]*p->vdir[i] -
631     duvw[2]*p->nrm[i] ;
632     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
633     return(0);
634     orig_dir += 2;
635     }
636     return(1);
637     }
638    
639     /* Reinhart/Tregenza sample generator */
640     int
641     sample_reinhart(PARAMS *p, int b, FVECT orig_dir[])
642     {
643     #define T_NALT 7
644     static const int tnaz[T_NALT] = {30, 30, 24, 24, 18, 12, 6};
645     const int RowMax = T_NALT*p->hsiz + 1;
646     const double RAH = (.5*PI)/(RowMax-.5);
647     #define rnaz(r) (r >= RowMax-1 ? 1 : p->hsiz*tnaz[r/p->hsiz])
648     int n = myRCmanager.accum;
649     int row, col;
650     double samp3[3];
651     double alt, azi;
652     double duvw[3];
653     int i;
654    
655     if (orig_dir == NULL) { /* just requesting number of bins? */
656     n = 0;
657     for (row = RowMax; row--; ) n += rnaz(row);
658     return(n);
659     }
660     row = 0; /* identify row & column */
661     col = b;
662     while (col >= rnaz(row)) {
663     col -= rnaz(row);
664     ++row;
665     }
666     while (n--) { /* stratified sampling */
667     SDmultiSamp(samp3, 3, (n+frandom())/myRCmanager.accum);
668     if (row >= RowMax-1) /* avoid crowding at zenith */
669     samp3[1] *= samp3[1];
670     alt = (row+samp3[1])*RAH;
671     azi = (2.*PI)*(col+samp3[2]-.5)/rnaz(row);
672     duvw[2] = cos(alt); /* measured from horizon */
673     duvw[0] = tsin(azi)*duvw[2];
674     duvw[1] = -tcos(azi)*duvw[2];
675     duvw[2] = sqrt(1. - duvw[2]*duvw[2]);
676     for (i = 3; i--; )
677     orig_dir[1][i] = -duvw[0]*p->udir[i] -
678     duvw[1]*p->vdir[i] -
679     duvw[2]*p->nrm[i] ;
680     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
681     return(0);
682     orig_dir += 2;
683     }
684     return(1);
685     #undef rnaz
686     #undef T_NALT
687     }
688    
689     /* Klems sample generator */
690     int
691     sample_klems(PARAMS *p, int b, FVECT orig_dir[])
692     {
693     static const char bname[4][20] = {
694     "LBNL/Klems Full",
695     "LBNL/Klems Half",
696     "INTERNAL ERROR",
697     "LBNL/Klems Quarter"
698     };
699     static ANGLE_BASIS *kbasis[4];
700     const int bi = p->hemis[1] - '1';
701     int n = myRCmanager.accum;
702     double samp2[2];
703     double duvw[3];
704     int i;
705    
706     if (!kbasis[bi]) { /* need to get basis, first */
707     for (i = 4; i--; )
708     if (!strcasecmp(abase_list[i].name, bname[bi])) {
709     kbasis[bi] = &abase_list[i];
710     break;
711     }
712     if (i < 0) {
713     sprintf(errmsg, "unknown hemisphere basis '%s'",
714     bname[bi]);
715     error(USER, errmsg);
716     return(0);
717     }
718     }
719     if (orig_dir == NULL) /* just requesting number of bins? */
720     return(kbasis[bi]->nangles);
721    
722     while (n--) { /* stratified sampling */
723     SDmultiSamp(samp2, 2, (n+frandom())/myRCmanager.accum);
724     if (!fo_getvec(duvw, b+samp2[1], kbasis[bi]))
725     return(0);
726     for (i = 3; i--; )
727     orig_dir[1][i] = -duvw[0]*p->udir[i] -
728     duvw[1]*p->vdir[i] -
729     duvw[2]*p->nrm[i] ;
730     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp2[0]))
731     return(0);
732     orig_dir += 2;
733     }
734     return(1);
735     }
736    
737     /* Prepare hemisphere basis sampler that will send rays to rcontrib */
738     int
739     prepare_sampler(PARAMS *p)
740     {
741     if (p->slist == NULL) { /* missing sample surface! */
742     error(USER, "no sender surface");
743     return(-1);
744     }
745     /* misplaced output file spec. */
746     if (p->outfn[0]) {
747     sprintf(errmsg, "ignoring output file in sender ('%s')",
748     p->outfn);
749     error(WARNING, errmsg);
750     }
751     /* check/set basis hemisphere */
752     if (!p->hemis[0]) {
753     error(USER, "missing sender sampling type");
754     return(-1);
755     }
756     if (normalize(p->nrm) == 0) {
757     error(USER, "undefined normal for sender sampling");
758     return(-1);
759     }
760     if (normalize(p->vup) == 0) {
761     if (fabs(p->nrm[2]) < .7)
762     p->vup[2] = 1;
763     else
764     p->vup[1] = 1;
765     }
766     fcross(p->udir, p->vup, p->nrm);
767     if (normalize(p->udir) == 0) {
768     error(USER, "up vector coincides with sender normal");
769     return(-1);
770     }
771     fcross(p->vdir, p->nrm, p->udir);
772     if (p->sign == '-') { /* left-handed coordinate system? */
773     p->udir[0] *= -1.;
774     p->udir[1] *= -1.;
775     p->udir[2] *= -1.;
776     }
777     if ((tolower(p->hemis[0]) == 'u') | (p->hemis[0] == '1'))
778     p->sample_basis = sample_uniform;
779     else if (tolower(p->hemis[0]) == 's' &&
780     tolower(p->hemis[1]) == 'c')
781     p->sample_basis = sample_shirchiu;
782     else if ((tolower(p->hemis[0]) == 'r') |
783     (tolower(p->hemis[0]) == 't'))
784     p->sample_basis = sample_reinhart;
785     else if (tolower(p->hemis[0]) == 'k') {
786     switch (p->hemis[1]) {
787     case '1':
788     case '2':
789     case '4':
790     break;
791     case 'f':
792     case 'F':
793     case '\0':
794     p->hemis[1] = '1';
795     break;
796     case 'h':
797     case 'H':
798     p->hemis[1] = '2';
799     break;
800     case 'q':
801     case 'Q':
802     p->hemis[1] = '4';
803     break;
804     default:
805     goto unrecognized;
806     }
807     p->hemis[2] = '\0';
808     p->sample_basis = sample_klems;
809     } else
810     goto unrecognized;
811     /* return number of bins */
812     return((*p->sample_basis)(p,0,NULL));
813     unrecognized:
814     sprintf(errmsg, "unrecognized sender sampling: h=%s", p->hemis);
815     error(USER, errmsg);
816     return(-1);
817     }
818    
819     /* Compute normal and area for polygon */
820     int
821     finish_polygon(SURF *p)
822     {
823     const int nv = p->nfargs / 3;
824     FVECT e1, e2, vc;
825     int i;
826    
827     memset(p->snrm, 0, sizeof(FVECT));
828     VSUB(e1, p->farg+3, p->farg);
829     for (i = 2; i < nv; i++) {
830     VSUB(e2, p->farg+3*i, p->farg);
831     VCROSS(vc, e1, e2);
832     p->snrm[0] += vc[0];
833     p->snrm[1] += vc[1];
834     p->snrm[2] += vc[2];
835     VCOPY(e1, e2);
836     }
837     p->area = normalize(p->snrm)*0.5;
838     return(p->area > FTINY*FTINY);
839     }
840    
841     /* Add a surface to our current parameters */
842     void
843     add_surface(int st, const char *oname, FILE *fp)
844     {
845     SURF *snew;
846     int n;
847     /* get floating-point arguments */
848     if (!fscanf(fp, "%d", &n)) return;
849     while (n-- > 0) fscanf(fp, "%*s");
850     if (!fscanf(fp, "%d", &n)) return;
851     while (n-- > 0) fscanf(fp, "%*d");
852     if (!fscanf(fp, "%d", &n) || n <= 0) return;
853     snew = (SURF *)malloc(sizeof(SURF) + sizeof(double)*(n-1));
854     if (snew == NULL)
855     error(SYSTEM, "out of memory in add_surface()");
856     strncpy(snew->sname, oname, sizeof(snew->sname)-1);
857     snew->sname[sizeof(snew->sname)-1] = '\0';
858     snew->styp = st;
859     snew->priv = NULL;
860     snew->nfargs = n;
861     for (n = 0; n < snew->nfargs; n++)
862     if (fscanf(fp, "%lf", &snew->farg[n]) != 1) {
863     sprintf(errmsg, "error reading arguments for '%s'",
864     oname);
865     error(USER, errmsg);
866     }
867     switch (st) {
868     case ST_RING:
869     if (snew->nfargs != 8)
870     goto badcount;
871     VCOPY(snew->snrm, snew->farg+3);
872     if (normalize(snew->snrm) == 0)
873     goto badnorm;
874     if (snew->farg[7] < snew->farg[6]) {
875     double t = snew->farg[7];
876     snew->farg[7] = snew->farg[6];
877     snew->farg[6] = t;
878     }
879     snew->area = PI*(snew->farg[7]*snew->farg[7] -
880     snew->farg[6]*snew->farg[6]);
881     break;
882     case ST_POLY:
883     if (snew->nfargs < 9 || snew->nfargs % 3)
884     goto badcount;
885     finish_polygon(snew);
886     break;
887     case ST_SOURCE:
888     if (snew->nfargs != 4)
889     goto badcount;
890     for (n = 3; n--; ) /* need to reverse "normal" */
891     snew->snrm[n] = -snew->farg[n];
892     if (normalize(snew->snrm) == 0)
893     goto badnorm;
894     snew->area = sin((PI/180./2.)*snew->farg[3]);
895     snew->area *= PI*snew->area;
896     break;
897     }
898     if (snew->area <= FTINY*FTINY) {
899     sprintf(errmsg, "zero area for surface '%s'", oname);
900     error(WARNING, errmsg);
901     free(snew);
902     return;
903     }
904     VSUM(curparams.nrm, curparams.nrm, snew->snrm, snew->area);
905     snew->next = curparams.slist;
906     curparams.slist = snew;
907     curparams.nsurfs++;
908     return;
909     badcount:
910     sprintf(errmsg, "bad argument count for surface element '%s'", oname);
911     error(USER, errmsg);
912     badnorm:
913     sprintf(errmsg, "bad orientation for surface element '%s'", oname);
914     error(USER, errmsg);
915     }
916    
917     /* Parse a receiver object (look for modifiers to add) */
918     int
919     add_recv_object(FILE *fp)
920     {
921     int st;
922     char thismod[128], otype[32], oname[128];
923     int n;
924    
925     if (fscanf(fp, "%s %s %s", thismod, otype, oname) != 3)
926     return(0); /* must have hit EOF! */
927     if (!strcmp(otype, "alias")) {
928     fscanf(fp, "%*s"); /* skip alias */
929     return(0);
930     }
931     /* is it a new receiver? */
932     if ((st = surf_type(otype)) != ST_NONE) {
933     if (strcmp(thismod, curmod)) {
934     if (curmod[0]) /* output last receiver? */
935     finish_receiver();
936     parse_params(&curparams, newparams);
937     newparams[0] = '\0';
938     strcpy(curmod, thismod);
939     }
940     add_surface(st, oname, fp); /* read & store surface */
941     return(1);
942     }
943     /* else skip arguments */
944     if (!fscanf(fp, "%d", &n)) return(0);
945     while (n-- > 0) fscanf(fp, "%*s");
946     if (!fscanf(fp, "%d", &n)) return(0);
947     while (n-- > 0) fscanf(fp, "%*d");
948     if (!fscanf(fp, "%d", &n)) return(0);
949     while (n-- > 0) fscanf(fp, "%*f");
950     return(0);
951     }
952    
953     /* Parse a sender object */
954     int
955     add_send_object(FILE *fp)
956     {
957     int st;
958     char thismod[128], otype[32], oname[128];
959     int n;
960    
961     if (fscanf(fp, "%s %s %s", thismod, otype, oname) != 3)
962     return(0); /* must have hit EOF! */
963     if (!strcmp(otype, "alias")) {
964     fscanf(fp, "%*s"); /* skip alias */
965     return(0);
966     }
967     /* is it a new surface? */
968     if ((st = surf_type(otype)) != ST_NONE) {
969     if (st == ST_SOURCE) {
970     error(USER, "cannot use source as a sender");
971     return(-1);
972     }
973     if (strcmp(thismod, curmod)) {
974     if (curmod[0])
975     error(WARNING, "multiple modifiers in sender");
976     strcpy(curmod, thismod);
977     }
978     parse_params(&curparams, newparams);
979     newparams[0] = '\0';
980     add_surface(st, oname, fp); /* read & store surface */
981     return(0);
982     }
983     /* else skip arguments */
984     if (!fscanf(fp, "%d", &n)) return(0);
985     while (n-- > 0) fscanf(fp, "%*s");
986     if (!fscanf(fp, "%d", &n)) return(0);
987     while (n-- > 0) fscanf(fp, "%*d");
988     if (!fscanf(fp, "%d", &n)) return(0);
989     while (n-- > 0) fscanf(fp, "%*f");
990     return(0);
991     }
992    
993     /* Load a Radiance scene using the given callback function for objects */
994     int
995     load_scene(const char *inspec, int (*ocb)(FILE *))
996     {
997     int rv = 0;
998     char inpbuf[1024];
999     FILE *fp;
1000     int c;
1001    
1002     if (*inspec == '!')
1003     fp = popen(inspec+1, "r");
1004     else
1005     fp = fopen(inspec, "r");
1006     if (fp == NULL) {
1007     sprintf(errmsg, "cannot load '%s'", inspec);
1008     error(SYSTEM, errmsg);
1009     return(-1);
1010     }
1011     while ((c = getc(fp)) != EOF) { /* load sender/receiver data */
1012     if (isspace(c)) /* skip leading white space */
1013     continue;
1014     if (c == '!') { /* read from a new command */
1015     inpbuf[0] = c;
1016     if (fgetline(inpbuf+1, sizeof(inpbuf)-1, fp) != NULL) {
1017     if ((c = load_scene(inpbuf, ocb)) < 0)
1018     return(c);
1019     rv += c;
1020     }
1021     continue;
1022     }
1023     if (c == '#') { /* parameters/comment */
1024     if ((c = getc(fp)) == EOF)
1025     break;
1026     if (c == '@' && fscanf(fp, "%s", inpbuf) == 1 &&
1027     (!strcmp(inpbuf, PARAMSNAME) ||
1028     !strcmp(inpbuf, progname))) {
1029     if (fgets(inpbuf, sizeof(inpbuf), fp) != NULL)
1030     strlcat(newparams, inpbuf, sizeof(newparams));
1031     continue;
1032     }
1033     while (c != '\n' && (c = getc(fp)) != EOF)
1034     ; /* ...skipping comment */
1035     continue;
1036     }
1037     ungetc(c, fp); /* else check object for receiver */
1038     c = (*ocb)(fp);
1039     if (c < 0)
1040     return(c);
1041     rv += c;
1042     }
1043     /* close our input stream */
1044     c = (*inspec == '!') ? pclose(fp) : fclose(fp);
1045     if (c != 0) {
1046     sprintf(errmsg, "error loading '%s'", inspec);
1047     error(SYSTEM, errmsg);
1048     return(-1);
1049     }
1050     return(rv);
1051     }
1052    
1053     void
1054     wputs( /* warning output function */
1055     const char *s
1056     )
1057     {
1058     if (erract[WARNING].pf == NULL)
1059     return;
1060     int lasterrno = errno;
1061     eputs(s);
1062     errno = lasterrno;
1063     }
1064    
1065     void
1066     eputs( /* put string to stderr */
1067     const char *s
1068     )
1069     {
1070     static int midline = 0;
1071    
1072     if (!*s)
1073     return;
1074     if (!midline++) {
1075     fputs(progname, stderr);
1076     fputs(": ", stderr);
1077     }
1078     fputs(s, stderr);
1079     if (s[strlen(s)-1] == '\n') {
1080     fflush(stderr);
1081     midline = 0;
1082     }
1083     }
1084    
1085     /* set input/output format */
1086     void
1087     setformat(const char *fmt)
1088     {
1089     switch (fmt[0]) {
1090     case 'f':
1091     case 'd':
1092     SET_FILE_BINARY(stdin);
1093     /* fall through */
1094     case 'a':
1095     inpfmt = fmt[0];
1096     break;
1097     case 'c':
1098     if (fmt[1])
1099     goto fmterr;
1100     outfmt = fmt[0]; // special case for output-only
1101     return;
1102     default:
1103     goto fmterr;
1104     }
1105     switch (fmt[1]) {
1106     case '\0':
1107     if (inpfmt == 'a')
1108     goto fmterr;
1109     outfmt = inpfmt;
1110     return;
1111     case 'f':
1112     case 'd':
1113     case 'c':
1114     outfmt = fmt[1];
1115     break;
1116     default:
1117     goto fmterr;
1118     }
1119     if (!fmt[2])
1120     return;
1121     fmterr:
1122     sprintf(errmsg, "unsupported i/o format: -f%s", fmt);
1123     error(USER, errmsg);
1124     }
1125    
1126     inline double
1127     pixjitter()
1128     {
1129     return(0.5 + dstrpix*(frandom()-0.5));
1130     }
1131    
1132     // Compute a set of view rays for the given pixel accumulator
1133     bool
1134     viewRays(FVECT orig_dir[], int x, int y)
1135     {
1136     for (int n = 0; n < myRCmanager.accum; orig_dir += 2, n++) {
1137     const double d = viewray(orig_dir[0], orig_dir[1], &ourview,
1138     (x+pixjitter())/myRCmanager.xres,
1139     (y+pixjitter())/myRCmanager.yres);
1140     // off-view sample?
1141     if (d < -FTINY || !jitteraperture(orig_dir[0], orig_dir[1],
1142     &ourview, dblur))
1143     memset(orig_dir+1, 0, sizeof(FVECT));
1144     else // else record distance to aft plane
1145     for (int i = 3*(d > FTINY); i--; )
1146     orig_dir[1][i] *= d;
1147     // accumulating identical samples?
1148     if (!n && (dstrpix <= 0.05) & (dblur <= FTINY)) {
1149     while (++n < myRCmanager.accum)
1150     memcpy(orig_dir[2*n], orig_dir[0], sizeof(FVECT)*2);
1151     break;
1152     }
1153     }
1154     return true;
1155     }
1156    
1157     // Load a set of rays for accumulation (do not normalize direction)
1158     int
1159     getRays(FVECT orig_dir[])
1160     {
1161     int n;
1162     // read directly if possible
1163     if (inpfmt == "_fd"[sizeof(RREAL)/sizeof(float)])
1164     return(getbinary(orig_dir[0], sizeof(FVECT)*2,
1165     myRCmanager.accum, stdin));
1166    
1167     for (n = 0; n < myRCmanager.accum; orig_dir += 2, n++) {
1168     switch (inpfmt) {
1169     #ifdef SMLFLT
1170     case 'd': { double dvin[6];
1171     if (getbinary(dvin, sizeof(dvin), 1, stdin) != 1)
1172     break;
1173     for (int i = 6; i--; ) orig_dir[0][i] = dvin[i];
1174     } continue;
1175     #else
1176     case 'f': { float fvin[6];
1177     if (getbinary(fvin, sizeof(fvin), 1, stdin) != 1)
1178     break;
1179     for (int i = 6; i--; ) orig_dir[0][i] = fvin[i];
1180     } continue;
1181     #endif
1182     case 'a':
1183     if (scanf(FVFORMAT, &orig_dir[0][0], &orig_dir[0][1],
1184     &orig_dir[0][2]) != 3)
1185     break;
1186     if (scanf(FVFORMAT, &orig_dir[1][0], &orig_dir[1][1],
1187     &orig_dir[1][2]) != 3)
1188     break;
1189     continue;
1190     }
1191     break;
1192     }
1193     return(n);
1194     }
1195    
1196     /* Set default options */
1197     void
1198     default_options()
1199     {
1200     rand_samp = 1;
1201     dstrsrc = 0.9;
1202     directrelay = 3;
1203     vspretest = 512;
1204     srcsizerat = .2;
1205     specthresh = .02;
1206     specjitter = 1.;
1207     maxdepth = -10;
1208     minweight = 2e-3;
1209     ambres = 256;
1210     ambdiv = 350;
1211     ambounce = 1;
1212     }
1213    
1214     /* Set overriding options */
1215     void
1216     override_options()
1217     {
1218     shadthresh = 0;
1219     ambssamp = 0;
1220     ambacc = 0;
1221     }
1222    
1223     void
1224     onsig( /* fatal signal */
1225     int signo
1226     )
1227     {
1228     static int gotsig = 0;
1229    
1230     if (gotsig++) /* two signals and we're gone! */
1231     _exit(signo);
1232    
1233     #ifdef SIGALRM
1234     alarm(600); /* allow 10 minutes to clean up */
1235     signal(SIGALRM, SIG_DFL); /* make certain we do die */
1236     #endif
1237     eputs("signal - ");
1238     eputs(sigerr[signo]);
1239     eputs("\n");
1240     quit(3);
1241     }
1242    
1243     void
1244     sigdie( /* set fatal signal */
1245     int signo,
1246     const char *msg
1247     )
1248     {
1249     if (signal(signo, onsig) == SIG_IGN)
1250     signal(signo, SIG_IGN);
1251     sigerr[signo] = msg;
1252     }
1253    
1254     /* Run rfluxmtx equivalent without leaning on r[x]contrib */
1255     int
1256     main(int argc, char *argv[])
1257     {
1258     #define check(ol,al) if (argv[a][ol] || \
1259     badarg(argc-a-1,argv+a+1,al)) \
1260     goto userr
1261     #define check_bool(olen,var) switch (argv[a][olen]) { \
1262     case '\0': var = !var; break; \
1263     case 'y': case 'Y': case 't': case 'T': \
1264     case '+': case '1': var = 1; break; \
1265     case 'n': case 'N': case 'f': case 'F': \
1266     case '-': case '0': var = 0; break; \
1267     default: goto userr; }
1268     bool gotView = false;
1269     double pixaspect = 1.;
1270     int nproc = 1;
1271     int nout = 0;
1272     char *outfn = NULL;
1273     FVECT *rayarr = NULL;
1274     const char *sendfn = "";
1275     PARAMS sendparams;
1276     int rval;
1277     int a, i;
1278     /* set global progname */
1279     fixargv0(argv[0]);
1280     // just asking version?
1281     if (argc == 2 && !strcmp(argv[1], "-version")) {
1282     puts(VersionID);
1283     quit(0);
1284     }
1285     initfunc(); /* initialize calcomp routines first */
1286     calcontext(RCCONTEXT);
1287     /* set rcontrib defaults */
1288     default_options();
1289     /* get command-line options */
1290     for (a = 1; a < argc-2; a++) {
1291     /* check for argument expansion */
1292     while ((rval = expandarg(&argc, &argv, a)) > 0)
1293     ;
1294     if (rval < 0) {
1295     sprintf(errmsg, "cannot expand '%s'", argv[a]);
1296     error(SYSTEM, errmsg);
1297     }
1298     if (argv[a][0] != '-' || !argv[a][1])
1299     break; // break from options
1300     rval = getrenderopt(argc-a, argv+a);
1301     if (rval >= 0) {
1302     a += rval;
1303     continue;
1304     }
1305     rval = getviewopt(&ourview, argc-a, argv+a);
1306     if (rval >= 0) {
1307     a += rval;
1308     gotView = true;
1309     continue;
1310     }
1311     switch (argv[a][1]) { /* !! Keep consistent !! */
1312     case 'W': /* verbose mode */
1313     check_bool(2,verby);
1314     break;
1315     case 'v': // view file
1316     if (argv[a][2] != 'f')
1317     goto userr;
1318     check(3,"s");
1319     rval = viewfile(argv[++a], &ourview, NULL);
1320     if (rval < 0) {
1321     sprintf(errmsg, "cannot open view file '%s'", argv[a]);
1322     error(SYSTEM, errmsg);
1323     } else if (!rval) {
1324     sprintf(errmsg, "bad view file '%s'", argv[a]);
1325     error(USER, errmsg);
1326     }
1327     gotView = true;
1328     break;
1329     case 'p': // -pj, -pd, or -pa
1330     switch (argv[a][2]) {
1331     case 'j': /* jitter */
1332     check(3,"f");
1333     dstrpix = atof(argv[++a]);
1334     break;
1335     case 'd': /* aperture */
1336     check(3,"f");
1337     dblur = atof(argv[++a]);
1338     break;
1339     case 'a': /* aspect */
1340     check(3,"f");
1341     pixaspect = atof(argv[++a]);
1342     break;
1343     default:
1344     goto userr;
1345     }
1346     break;
1347     case 'n': /* number of processes */
1348     check(2,"i");
1349     nproc = atoi(argv[++a]);
1350     if (nproc < 0 && (nproc += RadSimulManager::GetNCores()) <= 0)
1351     nproc = 1;
1352     break;
1353     case 'V': /* output contributions? */
1354     rval = myRCmanager.HasFlag(RCcontrib);
1355     check_bool(2,rval);
1356     myRCmanager.SetFlag(RCcontrib, rval);
1357     break;
1358     case 'x': /* x resolution */
1359     check(2,"i");
1360     myRCmanager.xres = atoi(argv[++a]);
1361     break;
1362     case 'y': /* y resolution */
1363     check(2,"i");
1364     myRCmanager.yres = atoi(argv[++a]);
1365     break;
1366     case 'w': /* warnings on/off */
1367     rval = (erract[WARNING].pf != NULL);
1368     check_bool(2,rval);
1369     if (rval) erract[WARNING].pf = wputs;
1370     else erract[WARNING].pf = NULL;
1371     break;
1372     case 'l': /* limit distance */
1373     if (argv[a][2] != 'd')
1374     goto userr;
1375     rval = myRCmanager.HasFlag(RTlimDist);
1376     check_bool(3,rval);
1377     myRCmanager.SetFlag(RTlimDist, rval);
1378     break;
1379     case 'I': /* immed. irradiance */
1380     rval = myRCmanager.HasFlag(RTimmIrrad);
1381     check_bool(2,rval);
1382     myRCmanager.SetFlag(RTimmIrrad, rval);
1383     break;
1384     case 'f': /* format */
1385     if (argv[a][2] == 'o')
1386     goto userr; /* -fo is not optional */
1387     setformat(argv[a]+2);
1388     myRCmanager.SetDataFormat(outfmt);
1389     break;
1390     case 'o': /* output file */
1391     check(2,"s");
1392     outfn = argv[++a];
1393     break;
1394     case 'c': /* sample count */
1395     check(2,"i");
1396     myRCmanager.accum = atoi(argv[++a]);
1397     break;
1398     default: /* anything else is verbotten */
1399     goto userr;
1400     }
1401     }
1402     if (a > argc-2)
1403     goto userr;
1404    
1405     override_options(); /* override critical options */
1406    
1407     if (!gotView)
1408     sendfn = argv[a++];
1409     else if (argv[a][0] == '-')
1410     error(USER, "view specification incompatible with pass-through mode");
1411     /* assign sender & receiver inputs */
1412     if (gotView) { // picture output?
1413     const char * err = setview(&ourview);
1414     if (err != NULL)
1415     error(USER, err);
1416     normaspect(viewaspect(&ourview), &pixaspect,
1417     &myRCmanager.xres, &myRCmanager.yres);
1418     if ((myRCmanager.xres <= 0) | (myRCmanager.yres <= 0))
1419     error(USER, "missing or illegal -x and -y resolution");
1420     myRCmanager.SetFlag(RTlimDist, ourview.vaft > FTINY);
1421     if (myRCmanager.accum <= 0)
1422     myRCmanager.accum = 1;
1423     } else if (sendfn[0] == '-') { // pass-through mode?
1424     if (sendfn[1]) goto userr;
1425     sendfn = NULL;
1426     if (myRCmanager.accum <= 0)
1427     myRCmanager.accum = 1;
1428     } else { // else surface sampling
1429     if (do_irrad | myRCmanager.HasFlag(RTimmIrrad))
1430     error(USER, "-i, -I not allowed in surface-sampling mode");
1431     myRCmanager.SetFlag(RTlimDist, false);
1432     if (load_scene(sendfn, add_send_object) < 0)
1433     quit(1);
1434     sendparams = curparams; // save sender params & clear current
1435     memset(&curparams, 0, sizeof(PARAMS));
1436     curmod[0] = '\0';
1437     myRCmanager.yres = prepare_sampler(&sendparams);
1438     if (myRCmanager.yres <= 0)
1439     quit(1);
1440     myRCmanager.xres = 0;
1441     if (myRCmanager.accum <= 1)
1442     myRCmanager.accum = 10000;
1443     }
1444     if (outfn != NULL) // saved from -o option above?
1445     strlcpy(curparams.outfn, outfn, sizeof(curparams.outfn));
1446     // get ready to rock...
1447     if (setspectrsamp(CNDX, WLPART) < 0)
1448     error(USER, "unsupported spectral sampling");
1449     /* set up signal handling */
1450     sigdie(SIGINT, "Interrupt");
1451     #ifdef SIGHUP
1452     sigdie(SIGHUP, "Hangup");
1453     #endif
1454     sigdie(SIGTERM, "Terminate");
1455     #ifdef SIGPIPE
1456     sigdie(SIGPIPE, "Broken pipe");
1457     #endif
1458     #ifdef SIGALRM
1459     sigdie(SIGALRM, "Alarm clock");
1460     #endif
1461     #ifdef SIGXCPU
1462     sigdie(SIGXCPU, "CPU limit exceeded");
1463     sigdie(SIGXFSZ, "File size exceeded");
1464     #endif
1465     #ifdef NICE
1466     nice(NICE); /* lower priority */
1467     #endif
1468     rayarr = new FVECT [myRCmanager.accum*2];
1469     // load octree
1470     if (!myRCmanager.LoadOctree(oconv_command(argc-a, argv+a)))
1471     quit(1);
1472     // add to header
1473     myRCmanager.AddHeader(argc, argv);
1474     {
1475     char buf[128] = "SOFTWARE= ";
1476     strcpy(buf+10, VersionID);
1477     myRCmanager.AddHeader(buf);
1478     if (gotView) {
1479     sprintf(buf, "%s%s", VIEWSTR, viewopt(&ourview));
1480     if (strlen(buf) > VIEWSTRL+3)
1481     myRCmanager.AddHeader(buf);
1482     }
1483     if ((pixaspect > .005) && (pixaspect < .995) |
1484     (pixaspect > 1.005)) {
1485     sprintf(buf, "%s%f", ASPECTSTR, pixaspect);
1486     myRCmanager.AddHeader(buf);
1487     }
1488     }
1489     /* assign receiver modifiers */
1490     if (load_scene(argv[a], add_recv_object) < 0)
1491     quit(1);
1492     finish_receiver(); // makes final AddModifier() call
1493     // prepare output(s)
1494     myRCmanager.outOp = RCOforce; // mandatory rcontrib -fo+ mode
1495     if (myRCmanager.PrepOutput() < 0)
1496     error(USER, "issue creating output file(s)");
1497     if (verby) // get output file count?
1498     for (const RcontribOutput *op = myRCmanager.GetOutput();
1499     op != NULL; op = op->Next())
1500     ++nout;
1501     if (nproc > 1) { // set #processes
1502     if (verby)
1503     fprintf(stderr, "%s: starting %d subprocesses\n", progname, nproc);
1504     myRCmanager.SetThreadCount(nproc);
1505     }
1506     if (gotView) { // picture generation mode?
1507     if (verby) {
1508     fprintf(stderr, "%s: computing %d %dx%d pictures",
1509     progname, nout, myRCmanager.xres, myRCmanager.yres);
1510     if (myRCmanager.accum > 1)
1511     fprintf(stderr, " with %d samples/pixel\n", myRCmanager.accum);
1512     else
1513     fputc('\n', stderr);
1514     }
1515     for (i = myRCmanager.yres; i--; ) // from the top!
1516     for (int x = 0; x < myRCmanager.xres; x++) {
1517     if (!viewRays(rayarr, x, i))
1518     quit(1);
1519     if (myRCmanager.ComputeRecord(rayarr) != myRCmanager.accum)
1520     error(USER, "failed call to ComputeRecord()");
1521     }
1522     } else if (sendfn == NULL) { // pass-through mode?
1523     #ifdef getc_unlocked
1524     flockfile(stdin);
1525     #endif
1526     if (verby)
1527     fprintf(stderr,
1528     "%s: computing %d rows in %d matrices with %d samples/row\n",
1529     progname, myRCmanager.GetRowMax(), nout, myRCmanager.accum);
1530     for (i = myRCmanager.GetRowMax(); i-- > 0; ) {
1531     if (getRays(rayarr) != myRCmanager.accum) {
1532     sprintf(errmsg, "ray read error after %d of %d",
1533     myRCmanager.GetRowCount(),
1534     myRCmanager.GetRowMax());
1535     error(USER, errmsg);
1536     }
1537     if (myRCmanager.ComputeRecord(rayarr) != myRCmanager.accum)
1538     error(USER, "failed call to ComputeRecord()");
1539     }
1540     } else { // else surface-sampling mode
1541     if (verby) {
1542     fprintf(stderr,
1543     "%s: sampling %d directions in %d matrices with %d samples/direction",
1544     progname, myRCmanager.yres, nout, myRCmanager.accum);
1545     if (sendparams.nsurfs > 1)
1546     fprintf(stderr, " (%d surface elements)\n", sendparams.nsurfs);
1547     else
1548     fputc('\n', stderr);
1549     }
1550     for (i = 0; i < myRCmanager.yres; i++) {
1551     if (!(*sendparams.sample_basis)(&sendparams, i, rayarr))
1552     quit(1);
1553     if (myRCmanager.ComputeRecord(rayarr) != myRCmanager.accum)
1554     error(USER, "failed call to ComputeRecord()");
1555     }
1556     clear_params(&sendparams);
1557     }
1558     delete [] rayarr;
1559     quit(0); /* flushes and waits on children */
1560     userr:
1561     if (a < argc-2)
1562     fprintf(stderr, "%s: unsupported option '%s'", progname, argv[a]);
1563     fprintf(stderr, "Usage: %s [-W] [rxcontrib options] { sender.rad | view | - } receiver.rad [-i system.oct] [system.rad ..]\n",
1564     progname);
1565     quit(1);
1566     }
1567    
1568     /* Exit program */
1569     void
1570     quit(
1571     int code
1572     )
1573     {
1574     myRCmanager.FlushQueue(); // leave nothing in queue
1575    
1576     exit(code);
1577     }