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

File Contents

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