ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/rt/rxfluxmtx.cpp
Revision: 2.7
Committed: Fri Oct 24 19:58:19 2025 UTC (43 hours, 32 minutes ago) by greg
Branch: MAIN
CVS Tags: HEAD
Changes since 2.6: +11 -8 lines
Log Message:
feat(rxfluxmtx): Minor tweaks to recover messages in -W+ mode

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rxfluxmtx.cpp,v 2.6 2025/10/24 18:37:05 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 <time.h>
16 #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 int report_intvl = 0; /* reporting interval (seconds) */
39
40 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 static PARAMS curparams;
86 static char curmod[MAXSTR];
87 static char newparams[1024];
88
89 typedef int SURFSAMP(FVECT, SURF *, double);
90
91 static SURFSAMP ssamp_bad, ssamp_poly, ssamp_ring;
92
93 static SURFSAMP *orig_in_surf[4] = {
94 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 static void
103 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 static int
123 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 static char *
136 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 static int
183 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 static int
223 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 static void
302 finish_receiver()
303 {
304 bool uniform = false;
305 const char *calfn = NULL;
306 char binv[64] = "0";
307 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 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 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 static void
416 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 static int
427 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 static int
436 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 static int
464 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 static int
477 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 static int
549 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 static int
590 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 static int
617 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 square2disk(duvw, (b/p->hsiz + samp3[1])/p->hsiz,
630 (b%p->hsiz + samp3[2])/p->hsiz);
631 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 static int
645 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 static int
695 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 static int
743 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 static int
825 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 static void
847 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 static int
923 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 static int
959 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 static int
999 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 static void
1091 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 static inline double
1131 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 viewRayBundle(FVECT orig_dir[], int x, int y)
1139 {
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 // 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 getRayBundle(FVECT orig_dir[])
1208 {
1209 // read directly if possible
1210 if (inpfmt == "_fd"[sizeof(RREAL)/sizeof(float)])
1211 return (getbinary(orig_dir[0], sizeof(FVECT)*2,
1212 myRCmanager.accum, stdin) == myRCmanager.accum);
1213
1214 for (int n = 0; n < myRCmanager.accum; orig_dir += 2, n++)
1215 switch (inpfmt) {
1216 #ifdef SMLFLT
1217 case 'd': { double dvin[6];
1218 if (getbinary(dvin, sizeof(dvin), 1, stdin) != 1)
1219 return false;
1220 for (int i = 6; i--; ) orig_dir[0][i] = dvin[i];
1221 } break;
1222 #else
1223 case 'f': { float fvin[6];
1224 if (getbinary(fvin, sizeof(fvin), 1, stdin) != 1)
1225 return false;
1226 for (int i = 6; i--; ) orig_dir[0][i] = fvin[i];
1227 } break;
1228 #endif
1229 case 'a':
1230 if (scanf(FVFORMAT, &orig_dir[0][0], &orig_dir[0][1],
1231 &orig_dir[0][2]) != 3)
1232 return false;
1233 if (scanf(FVFORMAT, &orig_dir[1][0], &orig_dir[1][1],
1234 &orig_dir[1][2]) != 3)
1235 return false;
1236 break;
1237 default:
1238 error(INTERNAL, "unsupported input format");
1239 }
1240
1241 return true;
1242 }
1243
1244 /* Set default options */
1245 static void
1246 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 static void
1264 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 // 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 /* 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 case '+': case '1': var = true; break; \
1335 case 'n': case 'N': case 'f': case 'F': \
1336 case '-': case '0': var = false; break; \
1337 default: goto userr; }
1338 bool force_open = false;
1339 bool recover = false;
1340 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 for (a = 1; a < argc; a++) {
1363 /* 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 case 'f': /* format or force overwrite */
1457 if (argv[a][2] == 'o') {
1458 check_bool(3,force_open);
1459 break;
1460 }
1461 setformat(argv[a]+2);
1462 break;
1463 case 'r': // recover flag
1464 check_bool(2,recover);
1465 break;
1466 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 case 't': /* reporting interval */
1475 check(2,"i");
1476 report_intvl = atoi(argv[++a]);
1477 break;
1478 default: /* anything else is verbotten */
1479 goto userr;
1480 }
1481 }
1482 if (a > argc-1)
1483 goto userr;
1484
1485 override_options(); /* override critical options */
1486
1487 if (!gotView) {
1488 sendfn = argv[a++];
1489 if (a > argc-1)
1490 goto userr;
1491 } else if (argv[a][0] == '-')
1492 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 if (myRCmanager.HasFlag(RTimmIrrad))
1499 error(USER, "-I+ not allowed in view generation mode");
1500 normaspect(viewaspect(&ourview), &pixaspect,
1501 &myRCmanager.xres, &myRCmanager.yres);
1502 if ((myRCmanager.xres <= 0) | (myRCmanager.yres <= 0))
1503 error(USER, "missing or illegal -x and -y resolution");
1504 myRCmanager.SetFlag(RTlimDist, ourview.vaft > FTINY);
1505 if (myRCmanager.accum <= 0)
1506 myRCmanager.accum = 1;
1507 } else if (sendfn[0] == '-') { // pass-through mode?
1508 if (sendfn[1]) goto userr;
1509 sendfn = NULL;
1510 if (myRCmanager.accum <= 0)
1511 myRCmanager.accum = 1;
1512 } else { // else surface sampling
1513 if (do_irrad | myRCmanager.HasFlag(RTimmIrrad))
1514 error(USER, "-i+, -I+ not allowed in surface-sampling mode");
1515 myRCmanager.SetFlag(RTlimDist, false);
1516 if (load_scene(sendfn, add_send_object) < 0)
1517 quit(1);
1518 sendparams = curparams; // save sender params & clear current
1519 memset(&curparams, 0, sizeof(PARAMS));
1520 curmod[0] = '\0';
1521 myRCmanager.yres = prepare_sampler(&sendparams);
1522 if (myRCmanager.yres <= 0)
1523 quit(1);
1524 myRCmanager.xres = 0;
1525 if (myRCmanager.accum <= 1)
1526 myRCmanager.accum = 10000;
1527 }
1528 if (outfn != NULL) // saved from -o option above?
1529 strlcpy(curparams.outfn, outfn, sizeof(curparams.outfn));
1530 // get ready to rock...
1531 if (setspectrsamp(CNDX, WLPART) < 0)
1532 error(USER, "unsupported spectral sampling");
1533 /* set up signal handling */
1534 sigdie(SIGINT, "Interrupt");
1535 #ifdef SIGHUP
1536 sigdie(SIGHUP, "Hangup");
1537 #endif
1538 sigdie(SIGTERM, "Terminate");
1539 #ifdef SIGPIPE
1540 sigdie(SIGPIPE, "Broken pipe");
1541 #endif
1542 #ifdef SIGALRM
1543 sigdie(SIGALRM, "Alarm clock");
1544 #endif
1545 #ifdef SIGXCPU
1546 sigdie(SIGXCPU, "CPU limit exceeded");
1547 sigdie(SIGXFSZ, "File size exceeded");
1548 #endif
1549 #ifdef NICE
1550 nice(NICE); /* lower priority */
1551 #endif
1552 rayarr = new FVECT [myRCmanager.accum*2];
1553 // load octree
1554 if (!myRCmanager.LoadOctree(oconv_command(argc-a, argv+a)))
1555 quit(1);
1556 // add to header
1557 myRCmanager.AddHeader(argc, argv);
1558 {
1559 char buf[128] = "SOFTWARE= ";
1560 strcpy(buf+10, VersionID);
1561 myRCmanager.AddHeader(buf);
1562 if (gotView) {
1563 sprintf(buf, "%s%s", VIEWSTR, viewopt(&ourview));
1564 if (strlen(buf) > VIEWSTRL+3)
1565 myRCmanager.AddHeader(buf);
1566 }
1567 if ((pixaspect > .005) & ((pixaspect < .995) |
1568 (pixaspect > 1.005))) {
1569 sprintf(buf, "%s%f", ASPECTSTR, pixaspect);
1570 myRCmanager.AddHeader(buf);
1571 }
1572 }
1573 // set output format
1574 myRCmanager.SetDataFormat(outfmt);
1575 /* assign receiver modifiers */
1576 if (load_scene(argv[a], add_recv_object) < 0)
1577 quit(1);
1578 finish_receiver(); // makes final AddModifier() call
1579 // prepare output files
1580 if (recover) {
1581 if (force_open) {
1582 error(WARNING, "-r+ mode overrides -fo+");
1583 force_open = false;
1584 }
1585 myRCmanager.outOp = RCOrecover;
1586 } else if (force_open)
1587 myRCmanager.outOp = RCOforce;
1588 else
1589 myRCmanager.outOp = RCOnew;
1590 // rval = # rows recovered
1591 rval = myRCmanager.PrepOutput();
1592 if (rval < 0) // PrepOutput() failure?
1593 error(USER, "issue loading or creating output");
1594 // in case output is complete
1595 if (rval >= myRCmanager.GetRowMax()) {
1596 error(WARNING, "nothing left to compute");
1597 quit(0);
1598 }
1599 if (verby) { // get output file count?
1600 if (rval > 0) {
1601 sprintf(errmsg, "recovered %d of %d rows\n",
1602 rval, myRCmanager.GetRowMax());
1603 eputs(errmsg);
1604 }
1605 for (const RcontribOutput *op = myRCmanager.GetOutput();
1606 op != NULL; op = op->Next())
1607 ++nout;
1608 }
1609 if (nproc > 1) { // set #processes
1610 if (verby) {
1611 sprintf(errmsg, "starting %d subprocesses\n", nproc);
1612 eputs(errmsg);
1613 }
1614 myRCmanager.SetThreadCount(nproc);
1615 }
1616 if (gotView) { // picture generation mode?
1617 i = myRCmanager.GetRowCount();
1618 if (verby) {
1619 sprintf(errmsg, "%s %d %dx%d pictures\n",
1620 i ? "completing" : "computing",
1621 nout, myRCmanager.xres, myRCmanager.yres);
1622 if (myRCmanager.accum > 1)
1623 sprintf(errmsg+strlen(errmsg)-1, " with %d samples/pixel\n",
1624 myRCmanager.accum);
1625 eputs(errmsg);
1626 }
1627 i /= myRCmanager.xres;
1628 int xstart = myRCmanager.GetRowCount() - i*myRCmanager.xres;
1629 i = myRCmanager.yres - i;
1630 while (i--) { // compute pixel rows from top down
1631 for (int x = xstart; x < myRCmanager.xres; x++) {
1632 report_progress();
1633 if (!viewRayBundle(rayarr, x, i))
1634 quit(1);
1635 if (myRCmanager.ComputeRecord(rayarr) != myRCmanager.accum)
1636 error(USER, "failed call to ComputeRecord()");
1637 }
1638 xstart = 0; // all rows after first start at x=0
1639 }
1640 } else if (sendfn == NULL) { // pass-through mode?
1641 #ifdef getc_unlocked
1642 flockfile(stdin);
1643 #endif
1644 // skip completed rows
1645 for (i = 0; i < myRCmanager.GetRowCount(); i++)
1646 if (!skipRayBundle()) {
1647 sprintf(errmsg, "read error from stdin at row %d", i);
1648 error(SYSTEM, errmsg);
1649 }
1650 if (verby) {
1651 sprintf(errmsg, "computing %d%s rows in %d matrices\n",
1652 myRCmanager.GetRowMax()-i, i ? " remaining" : "", nout);
1653 if (myRCmanager.accum > 1)
1654 sprintf(errmsg+strlen(errmsg)-1, " with %d samples/row\n",
1655 myRCmanager.accum);
1656 eputs(errmsg);
1657 }
1658 for ( ; i < myRCmanager.GetRowMax(); i++) {
1659 report_progress();
1660 if (!getRayBundle(rayarr)) {
1661 sprintf(errmsg, "read error from stdin at row %d of %d",
1662 i, myRCmanager.GetRowMax());
1663 error(USER, errmsg);
1664 }
1665 if (myRCmanager.ComputeRecord(rayarr) != myRCmanager.accum)
1666 error(USER, "failed call to ComputeRecord()");
1667 }
1668 } else { // else surface-sampling mode
1669 i = myRCmanager.GetRowCount();
1670 if (verby) {
1671 sprintf(errmsg, "sampling %d%s directions in %d matrices with %d samples/direction\n",
1672 myRCmanager.yres-i, i ? " remaining" : "",
1673 nout, myRCmanager.accum);
1674 if (sendparams.nsurfs > 1)
1675 sprintf(errmsg+strlen(errmsg)-1, " (%d surface elements)\n", sendparams.nsurfs);
1676 eputs(errmsg);
1677 }
1678 for ( ; i < myRCmanager.yres; i++) {
1679 report_progress();
1680 if (!(*sendparams.sample_basis)(&sendparams, i, rayarr))
1681 quit(1);
1682 if (myRCmanager.ComputeRecord(rayarr) != myRCmanager.accum)
1683 error(USER, "failed call to ComputeRecord()");
1684 }
1685 clear_params(&sendparams);
1686 }
1687 delete [] rayarr;
1688 myRCmanager.FlushQueue();
1689 report_progress((report_intvl > 0) | verby);
1690 quit(0); /* waits on any children */
1691 userr:
1692 if (a < argc && argv[a][0] == '-')
1693 fprintf(stderr, "%s: unsupported/misplaced option '%s'\n", progname, argv[a]);
1694 fprintf(stderr, "Usage: %s [-W] [rcontrib options] { sender.rad | view | - } receiver.rad [-i system.oct] [system.rad ..]\n",
1695 progname);
1696 quit(1);
1697 }
1698
1699 /* Exit program */
1700 void
1701 quit(
1702 int code
1703 )
1704 {
1705 if (!code)
1706 myRCmanager.ClearModifiers();
1707
1708 exit(code);
1709 }