ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/rt/rxfluxmtx.cpp
Revision: 2.5
Committed: Thu Oct 23 23:34:34 2025 UTC (12 days, 23 hours ago) by greg
Branch: MAIN
Changes since 2.4: +4 -3 lines
Log Message:
fix(rxfluxmtx): Fix seg fault with no arguments

File Contents

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