ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.30
Committed: Wed Jul 22 04:23:27 2015 UTC (8 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.29: +35 -24 lines
Log Message:
Improved matrix multiplication accuracy by forcing promotions to double

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.30 static const char RCSid[] = "$Id: rfluxmtx.c,v 2.29 2015/05/21 05:54:54 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Calculate flux transfer matrix or matrices using rcontrib
6     */
7    
8     #include "copyright.h"
9    
10     #include <ctype.h>
11     #include <stdlib.h>
12     #include "rtio.h"
13     #include "rtmath.h"
14 greg 2.12 #include "rtprocess.h"
15 greg 2.1 #include "bsdf.h"
16     #include "bsdf_m.h"
17     #include "random.h"
18     #include "triangulate.h"
19     #include "platform.h"
20    
21     #ifdef getc_unlocked /* avoid horrendous overhead of flockfile */
22     #undef getc
23     #define getc getc_unlocked
24     #endif
25    
26     #ifdef _WIN32
27 greg 2.2 #define SPECIALS " \t\"$*?"
28 greg 2.1 #define QUOTCHAR '"'
29     #else
30 greg 2.2 #define SPECIALS " \t\n'\"()${}*?[];|&"
31 greg 2.1 #define QUOTCHAR '\''
32     #define ALTQUOT '"'
33     #endif
34    
35     #define MAXRCARG 512
36    
37     char *progname; /* global argv[0] */
38    
39 greg 2.16 int verbose = 0; /* verbose mode (< 0 no warnings) */
40 greg 2.1
41     char *rcarg[MAXRCARG+1] = {"rcontrib", "-fo+"};
42     int nrcargs = 2;
43    
44     const char overflowerr[] = "%s: too many arguments!\n";
45    
46     #define CHECKARGC(n) if (nrcargs >= MAXRCARG-(n)) \
47     { fprintf(stderr, overflowerr, progname); exit(1); }
48    
49 greg 2.2 int sampcnt = 0; /* sample count (0==unset) */
50 greg 2.1
51 greg 2.2 char *reinhfn = "reinhartb.cal";
52     char *shirchiufn = "disk2square.cal";
53 greg 2.1 char *kfullfn = "klems_full.cal";
54     char *khalffn = "klems_half.cal";
55     char *kquarterfn = "klems_quarter.cal";
56    
57 greg 2.5 /* string indicating parameters */
58     const char PARAMSTART[] = "@rfluxmtx";
59 greg 2.1
60     /* surface type IDs */
61     #define ST_NONE 0
62     #define ST_POLY 1
63     #define ST_RING 2
64     #define ST_SOURCE 3
65    
66     typedef struct surf_s {
67     struct surf_s *next; /* next surface in list */
68     void *priv; /* private data (malloc'ed) */
69     char sname[32]; /* surface name */
70     FVECT snrm; /* surface normal */
71 greg 2.4 double area; /* surface area / proj. solid angle */
72 greg 2.1 short styp; /* surface type */
73 greg 2.2 short nfargs; /* number of real arguments */
74     double farg[1]; /* real values (extends struct) */
75 greg 2.1 } SURF; /* surface structure */
76    
77     typedef struct {
78     FVECT uva[2]; /* tangent axes */
79     int ntris; /* number of triangles */
80     struct ptri {
81     float afrac; /* fraction of total area */
82     short vndx[3]; /* vertex indices */
83     } tri[1]; /* triangle array (extends struct) */
84     } POLYTRIS; /* triangulated polygon */
85    
86     typedef struct param_s {
87 greg 2.26 char sign; /* '-' for axis reversal */
88     char hemis[31]; /* hemispherical sampling spec. */
89 greg 2.1 int hsiz; /* hemisphere basis size */
90     int nsurfs; /* number of surfaces */
91     SURF *slist; /* list of surfaces */
92     FVECT vup; /* up vector (zero if unset) */
93     FVECT nrm; /* average normal direction */
94 greg 2.26 FVECT udir, vdir; /* tangent axes */
95 greg 2.1 char *outfn; /* output file name (receiver) */
96     int (*sample_basis)(struct param_s *p, int, FILE *);
97     } PARAMS; /* sender/receiver parameters */
98    
99     PARAMS curparams;
100     char curmod[128];
101 greg 2.6 char newparams[1024];
102 greg 2.1
103     typedef int SURFSAMP(FVECT, SURF *, double);
104    
105     static SURFSAMP ssamp_bad, ssamp_poly, ssamp_ring;
106    
107     SURFSAMP *orig_in_surf[4] = {
108     ssamp_bad, ssamp_poly, ssamp_ring, ssamp_bad
109     };
110    
111     /* Clear parameter set */
112     static void
113     clear_params(PARAMS *p, int reset_only)
114     {
115     while (p->slist != NULL) {
116     SURF *sdel = p->slist;
117     p->slist = sdel->next;
118     if (sdel->priv != NULL)
119     free(sdel->priv);
120     free(sdel);
121     }
122     if (reset_only) {
123     p->nsurfs = 0;
124     memset(p->nrm, 0, sizeof(FVECT));
125     memset(p->vup, 0, sizeof(FVECT));
126     p->outfn = NULL;
127     return;
128     }
129 greg 2.6 memset(p, 0, sizeof(PARAMS));
130 greg 2.1 }
131    
132     /* Get surface type from name */
133     static int
134     surf_type(const char *otype)
135     {
136     if (!strcmp(otype, "polygon"))
137     return(ST_POLY);
138     if (!strcmp(otype, "ring"))
139     return(ST_RING);
140     if (!strcmp(otype, "source"))
141     return(ST_SOURCE);
142     return(ST_NONE);
143     }
144    
145 greg 2.30 /* Check if any of the characters in str2 are found in str1 */
146     static int
147     matchany(const char *str1, const char *str2)
148     {
149     while (*str1) {
150     const char *cp = str2;
151     while (*cp)
152     if (*cp++ == *str1)
153     return(*str1);
154     ++str1;
155     }
156     return(0);
157     }
158    
159 greg 2.1 /* Add arguments to oconv command */
160     static char *
161     oconv_command(int ac, char *av[])
162     {
163 greg 2.13 static char oconvbuf[2048] = "!oconv -f ";
164     char *cp = oconvbuf + 10;
165     char *recv = *av++;
166    
167     if (ac-- <= 0)
168     return(NULL);
169 greg 2.1 while (ac-- > 0) {
170 greg 2.13 strcpy(cp, *av++);
171     while (*cp) cp++;
172     *cp++ = ' ';
173 greg 2.30 if (cp >= oconvbuf+(sizeof(oconvbuf)-32))
174     goto overrun;
175     }
176     /* receiver goes last */
177     if (matchany(recv, SPECIALS)) {
178     *cp++ = QUOTCHAR;
179     while (*recv) {
180     if (cp >= oconvbuf+(sizeof(oconvbuf)-3))
181     goto overrun;
182     *cp++ = *recv++;
183 greg 2.1 }
184 greg 2.30 *cp++ = QUOTCHAR;
185     *cp = '\0';
186     } else
187     strcpy(cp, recv);
188 greg 2.1 return(oconvbuf);
189 greg 2.30 overrun:
190     fputs(progname, stderr);
191     fputs(": too many file arguments!\n", stderr);
192     exit(1);
193 greg 2.1 }
194    
195     /* Convert a set of arguments into a command line for pipe() or system() */
196     static char *
197     convert_commandline(char *cmd, const int len, char *av[])
198     {
199     int match;
200     char *cp;
201    
202     for (cp = cmd; *av != NULL; av++) {
203     const int n = strlen(*av);
204     if (cp+n >= cmd+(len-3)) {
205     fputs(progname, stderr);
206     return(NULL);
207     }
208 greg 2.30 if (matchany(*av, SPECIALS)) {
209 greg 2.1 const int quote =
210     #ifdef ALTQUOT
211 greg 2.30 strchr(*av, QUOTCHAR) ? ALTQUOT :
212 greg 2.1 #endif
213     QUOTCHAR;
214     *cp++ = quote;
215     strcpy(cp, *av);
216     cp += n;
217     *cp++ = quote;
218     } else {
219     strcpy(cp, *av);
220     cp += n;
221     }
222     *cp++ = ' ';
223     }
224     if (cp <= cmd)
225     return(NULL);
226     *--cp = '\0';
227     return(cmd);
228     }
229    
230     /* Open a pipe to/from a command given as an argument list */
231     static FILE *
232     popen_arglist(char *av[], char *mode)
233     {
234     char cmd[10240];
235    
236     if (!convert_commandline(cmd, sizeof(cmd), av)) {
237     fputs(progname, stderr);
238     fputs(": command line too long in popen_arglist()\n", stderr);
239     return(NULL);
240     }
241 greg 2.16 if (verbose > 0)
242 greg 2.1 fprintf(stderr, "%s: opening pipe %s: %s\n",
243     progname, (*mode=='w') ? "to" : "from", cmd);
244     return(popen(cmd, mode));
245     }
246    
247     #ifdef _WIN32
248     /* Execute system command (Windows version) */
249     static int
250     my_exec(char *av[])
251     {
252     char cmd[10240];
253    
254     if (!convert_commandline(cmd, sizeof(cmd), av)) {
255     fputs(progname, stderr);
256     fputs(": command line too long in my_exec()\n", stderr);
257     return(1);
258     }
259 greg 2.16 if (verbose > 0)
260 greg 2.1 fprintf(stderr, "%s: running: %s\n", progname, cmd);
261     return(system(cmd));
262     }
263     #else
264     /* Execute system command in our stead (Unix version) */
265     static int
266     my_exec(char *av[])
267     {
268     char *compath;
269    
270     if ((compath = getpath((char *)av[0], getenv("PATH"), X_OK)) == NULL) {
271     fprintf(stderr, "%s: cannot locate %s\n", progname, av[0]);
272     return(1);
273     }
274 greg 2.16 if (verbose > 0) {
275 greg 2.1 char cmd[4096];
276     if (!convert_commandline(cmd, sizeof(cmd), av))
277     strcpy(cmd, "COMMAND TOO LONG TO SHOW");
278     fprintf(stderr, "%s: running: %s\n", progname, cmd);
279     }
280     execv(compath, av); /* successful call never returns */
281     perror(compath);
282     return(1);
283     }
284     #endif
285    
286     /* Get normalized direction vector from string specification */
287     static int
288     get_direction(FVECT dv, const char *s)
289     {
290     int sign = 1;
291     int axis = 0;
292    
293     memset(dv, 0, sizeof(FVECT));
294     nextchar:
295     switch (*s) {
296     case '+':
297     ++s;
298     goto nextchar;
299     case '-':
300     sign = -sign;
301     ++s;
302     goto nextchar;
303     case 'z':
304     case 'Z':
305     ++axis;
306     case 'y':
307     case 'Y':
308     ++axis;
309     case 'x':
310     case 'X':
311     dv[axis] = sign;
312     return(!s[1] | isspace(s[1]));
313     default:
314     break;
315     }
316     #ifdef SMLFLT
317     if (sscanf(s, "%f,%f,%f", &dv[0], &dv[1], &dv[2]) != 3)
318     #else
319     if (sscanf(s, "%lf,%lf,%lf", &dv[0], &dv[1], &dv[2]) != 3)
320     #endif
321     return(0);
322     dv[0] *= (RREAL)sign;
323     return(normalize(dv) > 0);
324     }
325    
326     /* Parse program parameters (directives) */
327     static int
328 greg 2.6 parse_params(PARAMS *p, char *pargs)
329 greg 2.1 {
330     char *cp = pargs;
331     int nparams = 0;
332     int i;
333    
334 greg 2.15 for ( ; ; ) {
335 greg 2.1 switch (*cp++) {
336     case 'h':
337     if (*cp++ != '=')
338     break;
339 greg 2.26 if ((*cp == '+') | (*cp == '-'))
340     p->sign = *cp++;
341     else
342     p->sign = '+';
343 greg 2.6 p->hsiz = 0;
344 greg 2.1 i = 0;
345     while (*cp && !isspace(*cp)) {
346     if (isdigit(*cp))
347 greg 2.6 p->hsiz = 10*p->hsiz + *cp - '0';
348     p->hemis[i++] = *cp++;
349 greg 2.1 }
350     if (!i)
351     break;
352 greg 2.6 p->hemis[i] = '\0';
353     p->hsiz += !p->hsiz;
354 greg 2.1 ++nparams;
355     continue;
356     case 'u':
357     if (*cp++ != '=')
358     break;
359 greg 2.6 if (!get_direction(p->vup, cp))
360 greg 2.1 break;
361 greg 2.15 while (*cp && !isspace(*cp++))
362     ;
363 greg 2.1 ++nparams;
364     continue;
365     case 'o':
366     if (*cp++ != '=')
367     break;
368     i = 0;
369     while (*cp && !isspace(*cp++))
370     i++;
371     if (!i)
372     break;
373     *--cp = '\0';
374 greg 2.6 p->outfn = savqstr(cp-i);
375 greg 2.1 *cp++ = ' ';
376     ++nparams;
377     continue;
378     case ' ':
379     case '\t':
380     case '\r':
381 greg 2.21 case '\n':
382 greg 2.15 continue;
383 greg 2.1 case '\0':
384     return(nparams);
385     default:
386     break;
387     }
388 greg 2.15 break;
389     }
390     fprintf(stderr, "%s: bad parameter string: %s", progname, pargs);
391 greg 2.1 exit(1);
392     return(-1); /* pro forma return */
393     }
394    
395     /* Add receiver arguments (directives) corresponding to the current modifier */
396     static void
397     finish_receiver(void)
398     {
399 greg 2.4 char sbuf[256];
400     int uniform = 0;
401 greg 2.2 char *calfn = NULL;
402 greg 2.3 char *params = NULL;
403 greg 2.2 char *binv = NULL;
404     char *binf = NULL;
405     char *nbins = NULL;
406 greg 2.1
407     if (!curmod[0]) {
408     fputs(progname, stderr);
409     fputs(": missing receiver surface!\n", stderr);
410     exit(1);
411     }
412     if (curparams.outfn != NULL) { /* add output file spec. */
413     CHECKARGC(2);
414     rcarg[nrcargs++] = "-o";
415     rcarg[nrcargs++] = curparams.outfn;
416     }
417 greg 2.4 /* check arguments */
418 greg 2.1 if (!curparams.hemis[0]) {
419     fputs(progname, stderr);
420     fputs(": missing hemisphere sampling type!\n", stderr);
421     exit(1);
422     }
423     if (normalize(curparams.nrm) == 0) {
424     fputs(progname, stderr);
425     fputs(": undefined normal for hemisphere sampling\n", stderr);
426     exit(1);
427     }
428 greg 2.10 if (normalize(curparams.vup) == 0) {
429 greg 2.1 if (fabs(curparams.nrm[2]) < .7)
430     curparams.vup[2] = 1;
431     else
432     curparams.vup[1] = 1;
433 greg 2.10 }
434 greg 2.4 /* determine sample type/bin */
435 greg 2.1 if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1') {
436 greg 2.20 sprintf(sbuf, "if(-Dx*%g-Dy*%g-Dz*%g,0,-1)",
437     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2]);
438     binv = savqstr(sbuf);
439     nbins = "1"; /* uniform sampling -- one bin */
440 greg 2.4 uniform = 1;
441 greg 2.1 } else if (tolower(curparams.hemis[0]) == 's' &&
442     tolower(curparams.hemis[1]) == 'c') {
443 greg 2.2 /* assign parameters */
444 greg 2.1 if (curparams.hsiz <= 1) {
445     fputs(progname, stderr);
446     fputs(": missing size for Shirley-Chiu sampling!\n", stderr);
447     exit(1);
448     }
449 greg 2.3 calfn = shirchiufn; shirchiufn = NULL;
450 greg 2.26 sprintf(sbuf, "SCdim=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g,RHS=%c1",
451 greg 2.2 curparams.hsiz,
452     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
453 greg 2.26 curparams.vup[0], curparams.vup[1], curparams.vup[2],
454     curparams.sign);
455 greg 2.3 params = savqstr(sbuf);
456 greg 2.2 binv = "scbin";
457     nbins = "SCdim*SCdim";
458 greg 2.1 } else if ((tolower(curparams.hemis[0]) == 'r') |
459     (tolower(curparams.hemis[0]) == 't')) {
460 greg 2.2 calfn = reinhfn; reinhfn = NULL;
461 greg 2.26 sprintf(sbuf, "MF=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g,RHS=%c1",
462 greg 2.3 curparams.hsiz,
463     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
464 greg 2.26 curparams.vup[0], curparams.vup[1], curparams.vup[2],
465     curparams.sign);
466 greg 2.3 params = savqstr(sbuf);
467     binv = "rbin";
468     nbins = "Nrbins";
469 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
470     !curparams.hemis[1] |
471     (tolower(curparams.hemis[1]) == 'f') |
472     (curparams.hemis[1] == '1')) {
473 greg 2.2 calfn = kfullfn; kfullfn = NULL;
474     binf = "kbin";
475     nbins = "Nkbins";
476 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
477     (tolower(curparams.hemis[1]) == 'h') |
478     (curparams.hemis[1] == '2')) {
479 greg 2.2 calfn = khalffn; khalffn = NULL;
480     binf = "khbin";
481     nbins = "Nkhbins";
482 greg 2.1 } else if (tolower(curparams.hemis[0]) == 'k' &&
483     (tolower(curparams.hemis[1]) == 'q') |
484     (curparams.hemis[1] == '4')) {
485 greg 2.2 calfn = kquarterfn; kquarterfn = NULL;
486     binf = "kqbin";
487     nbins = "Nkqbins";
488 greg 2.1 } else {
489     fprintf(stderr, "%s: unrecognized hemisphere sampling: h=%s\n",
490     progname, curparams.hemis);
491     exit(1);
492     }
493 greg 2.26 if (tolower(curparams.hemis[0]) == 'k') {
494     sprintf(sbuf, "RHS=%c1", curparams.sign);
495     params = savqstr(sbuf);
496     }
497 greg 2.4 if (!uniform & (curparams.slist->styp == ST_SOURCE)) {
498     SURF *sp;
499     for (sp = curparams.slist; sp != NULL; sp = sp->next)
500     if (fabs(sp->area - PI) > 1e-3) {
501     fprintf(stderr, "%s: source '%s' must be 180-degrees\n",
502     progname, sp->sname);
503     exit(1);
504     }
505     }
506 greg 2.2 if (calfn != NULL) { /* add cal file if needed */
507     CHECKARGC(2);
508     rcarg[nrcargs++] = "-f";
509     rcarg[nrcargs++] = calfn;
510     }
511 greg 2.3 if (params != NULL) { /* parameters _after_ cal file */
512     CHECKARGC(2);
513     rcarg[nrcargs++] = "-p";
514     rcarg[nrcargs++] = params;
515     }
516 greg 2.2 if (nbins != NULL) { /* add #bins if set */
517     CHECKARGC(2);
518     rcarg[nrcargs++] = "-bn";
519     rcarg[nrcargs++] = nbins;
520     }
521 greg 2.3 if (binv != NULL) {
522 greg 2.2 CHECKARGC(2); /* assign bin variable */
523     rcarg[nrcargs++] = "-b";
524     rcarg[nrcargs++] = binv;
525     } else if (binf != NULL) {
526     CHECKARGC(2); /* assign bin function */
527 greg 2.3 rcarg[nrcargs++] = "-b";
528 greg 2.2 sprintf(sbuf, "%s(%g,%g,%g,%g,%g,%g)", binf,
529     curparams.nrm[0], curparams.nrm[1], curparams.nrm[2],
530     curparams.vup[0], curparams.vup[1], curparams.vup[2]);
531     rcarg[nrcargs++] = savqstr(sbuf);
532     }
533     CHECKARGC(2); /* modifier argument goes last */
534 greg 2.1 rcarg[nrcargs++] = "-m";
535     rcarg[nrcargs++] = savqstr(curmod);
536     }
537    
538     /* Make randomly oriented tangent plane axes for given normal direction */
539     static void
540     make_axes(FVECT uva[2], const FVECT nrm)
541     {
542     int i;
543    
544 greg 2.29 if (!getperpendicular(uva[0], nrm, 1)) {
545 greg 2.1 fputs(progname, stderr);
546     fputs(": bad surface normal in make_axes!\n", stderr);
547     exit(1);
548     }
549 greg 2.19 fcross(uva[1], nrm, uva[0]);
550 greg 2.1 }
551    
552     /* Illegal sender surfaces end up here */
553     static int
554     ssamp_bad(FVECT orig, SURF *sp, double x)
555     {
556     fprintf(stderr, "%s: illegal sender surface '%s'\n",
557     progname, sp->sname);
558     return(0);
559     }
560    
561     /* Generate origin on ring surface from uniform random variable */
562     static int
563     ssamp_ring(FVECT orig, SURF *sp, double x)
564     {
565     FVECT *uva = (FVECT *)sp->priv;
566     double samp2[2];
567     double uv[2];
568     int i;
569    
570     if (uva == NULL) { /* need tangent axes */
571     uva = (FVECT *)malloc(sizeof(FVECT)*2);
572     if (uva == NULL) {
573     fputs(progname, stderr);
574     fputs(": out of memory in ssamp_ring!\n", stderr);
575     return(0);
576     }
577     make_axes(uva, sp->snrm);
578     sp->priv = (void *)uva;
579     }
580     SDmultiSamp(samp2, 2, x);
581 greg 2.6 samp2[0] = sqrt(samp2[0]*sp->area*(1./PI) + sp->farg[6]*sp->farg[6]);
582 greg 2.1 samp2[1] *= 2.*PI;
583     uv[0] = samp2[0]*tcos(samp2[1]);
584     uv[1] = samp2[0]*tsin(samp2[1]);
585     for (i = 3; i--; )
586     orig[i] = sp->farg[i] + uv[0]*uva[0][i] + uv[1]*uva[1][i];
587     return(1);
588     }
589    
590     /* Add triangle to polygon's list (call-back function) */
591     static int
592     add_triangle(const Vert2_list *tp, int a, int b, int c)
593     {
594     POLYTRIS *ptp = (POLYTRIS *)tp->p;
595     struct ptri *trip = ptp->tri + ptp->ntris++;
596    
597     trip->vndx[0] = a;
598     trip->vndx[1] = b;
599     trip->vndx[2] = c;
600     return(1);
601     }
602    
603     /* Generate origin on polygon surface from uniform random variable */
604     static int
605     ssamp_poly(FVECT orig, SURF *sp, double x)
606     {
607     POLYTRIS *ptp = (POLYTRIS *)sp->priv;
608     double samp2[2];
609     double *v0, *v1, *v2;
610     int i;
611    
612     if (ptp == NULL) { /* need to triangulate */
613     ptp = (POLYTRIS *)malloc(sizeof(POLYTRIS) +
614     sizeof(struct ptri)*(sp->nfargs/3 - 3));
615     if (ptp == NULL)
616     goto memerr;
617     if (sp->nfargs == 3) { /* simple case */
618     ptp->ntris = 1;
619     ptp->tri[0].vndx[0] = 0;
620     ptp->tri[0].vndx[1] = 1;
621     ptp->tri[0].vndx[2] = 2;
622     ptp->tri[0].afrac = 1;
623     } else {
624     Vert2_list *v2l = polyAlloc(sp->nfargs/3);
625     if (v2l == NULL)
626     goto memerr;
627     make_axes(ptp->uva, sp->snrm);
628     for (i = v2l->nv; i--; ) {
629     v2l->v[i].mX = DOT(sp->farg+3*i, ptp->uva[0]);
630     v2l->v[i].mY = DOT(sp->farg+3*i, ptp->uva[1]);
631     }
632     ptp->ntris = 0;
633     v2l->p = (void *)ptp;
634 greg 2.7 if (!polyTriangulate(v2l, add_triangle)) {
635     fprintf(stderr,
636     "%s: cannot triangulate polygon '%s'\n",
637     progname, sp->sname);
638 greg 2.1 return(0);
639 greg 2.7 }
640 greg 2.1 for (i = ptp->ntris; i--; ) {
641     int a = ptp->tri[i].vndx[0];
642     int b = ptp->tri[i].vndx[1];
643     int c = ptp->tri[i].vndx[2];
644     ptp->tri[i].afrac =
645     (v2l->v[a].mX*v2l->v[b].mY -
646     v2l->v[b].mX*v2l->v[a].mY +
647     v2l->v[b].mX*v2l->v[c].mY -
648     v2l->v[c].mX*v2l->v[b].mY +
649     v2l->v[c].mX*v2l->v[a].mY -
650     v2l->v[a].mX*v2l->v[c].mY) /
651     (2.*sp->area);
652     }
653     polyFree(v2l);
654     }
655     sp->priv = (void *)ptp;
656     }
657     /* pick triangle by partial area */
658     for (i = 0; i < ptp->ntris && x > ptp->tri[i].afrac; i++)
659     x -= ptp->tri[i].afrac;
660     SDmultiSamp(samp2, 2, x/ptp->tri[i].afrac);
661     samp2[0] *= samp2[1] = sqrt(samp2[1]);
662     samp2[1] = 1. - samp2[1];
663     v0 = sp->farg + 3*ptp->tri[i].vndx[0];
664     v1 = sp->farg + 3*ptp->tri[i].vndx[1];
665     v2 = sp->farg + 3*ptp->tri[i].vndx[2];
666     for (i = 3; i--; )
667     orig[i] = v0[i] + samp2[0]*(v1[i] - v0[i])
668     + samp2[1]*(v2[i] - v0[i]) ;
669     return(1);
670     memerr:
671     fputs(progname, stderr);
672     fputs(": out of memory in ssamp_poly!\n", stderr);
673     return(0);
674     }
675    
676     /* Compute sample origin based on projected areas of sender subsurfaces */
677     static int
678     sample_origin(PARAMS *p, FVECT orig, const FVECT rdir, double x)
679     {
680     static double *projsa;
681     static int nall;
682     double tarea = 0;
683     int i;
684     SURF *sp;
685     /* special case for lone surface */
686     if (p->nsurfs == 1) {
687     sp = p->slist;
688 greg 2.7 if (DOT(sp->snrm, rdir) >= -FTINY) {
689     fprintf(stderr,
690     "%s: internal - sample behind sender '%s'\n",
691     progname, sp->sname);
692     return(0);
693     }
694 greg 2.1 return((*orig_in_surf[sp->styp])(orig, sp, x));
695     }
696     if (p->nsurfs > nall) { /* (re)allocate surface area cache */
697     if (projsa) free(projsa);
698     projsa = (double *)malloc(sizeof(double)*p->nsurfs);
699     if (!projsa) return(0);
700     nall = p->nsurfs;
701     }
702     /* compute projected areas */
703     for (i = 0, sp = p->slist; sp != NULL; i++, sp = sp->next) {
704     projsa[i] = -DOT(sp->snrm, rdir) * sp->area;
705     tarea += projsa[i] *= (double)(projsa[i] > FTINY);
706     }
707 greg 2.7 if (tarea <= FTINY) { /* wrong side of sender? */
708     fputs(progname, stderr);
709     fputs(": internal - sample behind all sender elements!\n",
710     stderr);
711 greg 2.1 return(0);
712 greg 2.7 }
713 greg 2.1 tarea *= x; /* get surface from list */
714     for (i = 0, sp = p->slist; tarea > projsa[i]; sp = sp->next)
715     tarea -= projsa[i++];
716     return((*orig_in_surf[sp->styp])(orig, sp, tarea/projsa[i]));
717     }
718    
719     /* Uniform sample generator */
720     static int
721     sample_uniform(PARAMS *p, int b, FILE *fp)
722     {
723     int n = sampcnt;
724     double samp3[3];
725     double duvw[3];
726     FVECT orig_dir[2];
727     int i;
728    
729     if (fp == NULL) /* just requesting number of bins? */
730     return(1);
731    
732     while (n--) { /* stratified hemisphere sampling */
733     SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
734     SDsquare2disk(duvw, samp3[1], samp3[2]);
735     duvw[2] = -sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
736     for (i = 3; i--; )
737     orig_dir[1][i] = duvw[0]*p->udir[i] +
738     duvw[1]*p->vdir[i] +
739     duvw[2]*p->nrm[i] ;
740     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
741     return(0);
742     if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
743     return(0);
744     }
745     return(1);
746     }
747    
748     /* Shirly-Chiu sample generator */
749     static int
750     sample_shirchiu(PARAMS *p, int b, FILE *fp)
751     {
752     int n = sampcnt;
753     double samp3[3];
754     double duvw[3];
755     FVECT orig_dir[2];
756     int i;
757    
758     if (fp == NULL) /* just requesting number of bins? */
759     return(p->hsiz*p->hsiz);
760    
761     while (n--) { /* stratified sampling */
762     SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
763     SDsquare2disk(duvw, (b/p->hsiz + samp3[1])/curparams.hsiz,
764     (b%p->hsiz + samp3[2])/curparams.hsiz);
765     duvw[2] = sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]);
766     for (i = 3; i--; )
767     orig_dir[1][i] = -duvw[0]*p->udir[i] -
768     duvw[1]*p->vdir[i] -
769     duvw[2]*p->nrm[i] ;
770     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
771     return(0);
772     if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
773     return(0);
774     }
775     return(1);
776     }
777    
778     /* Reinhart/Tregenza sample generator */
779     static int
780     sample_reinhart(PARAMS *p, int b, FILE *fp)
781     {
782     #define T_NALT 7
783     static const int tnaz[T_NALT] = {30, 30, 24, 24, 18, 12, 6};
784     const int RowMax = T_NALT*p->hsiz + 1;
785 greg 2.7 const double RAH = (.5*PI)/(RowMax-.5);
786 greg 2.1 #define rnaz(r) (r >= RowMax-1 ? 1 : p->hsiz*tnaz[r/p->hsiz])
787     int n = sampcnt;
788     int row, col;
789     double samp3[3];
790     double alt, azi;
791     double duvw[3];
792     FVECT orig_dir[2];
793     int i;
794    
795     if (fp == NULL) { /* just requesting number of bins? */
796     n = 0;
797     for (row = RowMax; row--; ) n += rnaz(row);
798     return(n);
799     }
800     row = 0; /* identify row & column */
801     col = b;
802     while (col >= rnaz(row)) {
803     col -= rnaz(row);
804     ++row;
805     }
806     while (n--) { /* stratified sampling */
807     SDmultiSamp(samp3, 3, (n+frandom())/sampcnt);
808     alt = (row+samp3[1])*RAH;
809     azi = (2.*PI)*(col+samp3[2]-.5)/rnaz(row);
810 greg 2.7 duvw[2] = cos(alt); /* measured from horizon */
811 greg 2.18 duvw[0] = tsin(azi)*duvw[2];
812     duvw[1] = tcos(azi)*duvw[2];
813 greg 2.1 duvw[2] = sqrt(1. - duvw[2]*duvw[2]);
814     for (i = 3; i--; )
815     orig_dir[1][i] = -duvw[0]*p->udir[i] -
816     duvw[1]*p->vdir[i] -
817     duvw[2]*p->nrm[i] ;
818     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0]))
819     return(0);
820     if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
821     return(0);
822     }
823     return(1);
824     #undef rnaz
825     #undef T_NALT
826     }
827    
828     /* Klems sample generator */
829     static int
830     sample_klems(PARAMS *p, int b, FILE *fp)
831     {
832     static const char bname[4][20] = {
833     "LBNL/Klems Full",
834     "LBNL/Klems Half",
835     "INTERNAL ERROR",
836     "LBNL/Klems Quarter"
837     };
838     static ANGLE_BASIS *kbasis[4];
839     const int bi = p->hemis[1] - '1';
840     int n = sampcnt;
841     double samp2[2];
842     double duvw[3];
843     FVECT orig_dir[2];
844     int i;
845    
846     if (!kbasis[bi]) { /* need to get basis, first */
847     for (i = 4; i--; )
848     if (!strcasecmp(abase_list[i].name, bname[bi])) {
849     kbasis[bi] = &abase_list[i];
850     break;
851     }
852     if (i < 0) {
853     fprintf(stderr, "%s: unknown hemisphere basis '%s'\n",
854     progname, bname[bi]);
855     return(0);
856     }
857     }
858     if (fp == NULL) /* just requesting number of bins? */
859     return(kbasis[bi]->nangles);
860    
861     while (n--) { /* stratified sampling */
862     SDmultiSamp(samp2, 2, (n+frandom())/sampcnt);
863 greg 2.27 if (!fo_getvec(duvw, b+samp2[1], kbasis[bi]))
864 greg 2.1 return(0);
865     for (i = 3; i--; )
866 greg 2.26 orig_dir[1][i] = -duvw[0]*p->udir[i] -
867     duvw[1]*p->vdir[i] -
868 greg 2.1 duvw[2]*p->nrm[i] ;
869     if (!sample_origin(p, orig_dir[0], orig_dir[1], samp2[0]))
870     return(0);
871     if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2)
872     return(0);
873     }
874     return(1);
875     }
876    
877     /* Prepare hemisphere basis sampler that will send rays to rcontrib */
878     static int
879     prepare_sampler(void)
880     {
881     if (curparams.slist == NULL) { /* missing sample surface! */
882     fputs(progname, stderr);
883     fputs(": no sender surface!\n", stderr);
884     return(-1);
885     }
886 greg 2.16 /* misplaced output file spec. */
887     if ((curparams.outfn != NULL) & (verbose >= 0))
888 greg 2.1 fprintf(stderr, "%s: warning - ignoring output file in sender ('%s')\n",
889     progname, curparams.outfn);
890     /* check/set basis hemisphere */
891     if (!curparams.hemis[0]) {
892     fputs(progname, stderr);
893     fputs(": missing sender sampling type!\n", stderr);
894     return(-1);
895     }
896     if (normalize(curparams.nrm) == 0) {
897     fputs(progname, stderr);
898     fputs(": undefined normal for sender sampling\n", stderr);
899     return(-1);
900     }
901 greg 2.10 if (normalize(curparams.vup) == 0) {
902 greg 2.1 if (fabs(curparams.nrm[2]) < .7)
903     curparams.vup[2] = 1;
904     else
905     curparams.vup[1] = 1;
906 greg 2.10 }
907 greg 2.26 fcross(curparams.udir, curparams.vup, curparams.nrm);
908 greg 2.1 if (normalize(curparams.udir) == 0) {
909     fputs(progname, stderr);
910     fputs(": up vector coincides with sender normal\n", stderr);
911     return(-1);
912     }
913 greg 2.26 fcross(curparams.vdir, curparams.nrm, curparams.udir);
914     if (curparams.sign == '-') { /* left-handed coordinate system? */
915     curparams.udir[0] *= -1.;
916     curparams.udir[1] *= -1.;
917     curparams.udir[2] *= -1.;
918     }
919 greg 2.1 if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1')
920     curparams.sample_basis = sample_uniform;
921     else if (tolower(curparams.hemis[0]) == 's' &&
922     tolower(curparams.hemis[1]) == 'c')
923     curparams.sample_basis = sample_shirchiu;
924     else if ((tolower(curparams.hemis[0]) == 'r') |
925     (tolower(curparams.hemis[0]) == 't'))
926     curparams.sample_basis = sample_reinhart;
927     else if (tolower(curparams.hemis[0]) == 'k') {
928     switch (curparams.hemis[1]) {
929     case '1':
930     case '2':
931     case '4':
932     break;
933     case 'f':
934     case 'F':
935     case '\0':
936     curparams.hemis[1] = '1';
937     break;
938     case 'h':
939     case 'H':
940     curparams.hemis[1] = '2';
941     break;
942     case 'q':
943     case 'Q':
944     curparams.hemis[1] = '4';
945     break;
946     default:
947     goto unrecognized;
948     }
949     curparams.hemis[2] = '\0';
950     curparams.sample_basis = sample_klems;
951     } else
952     goto unrecognized;
953     /* return number of bins */
954     return((*curparams.sample_basis)(&curparams,0,NULL));
955     unrecognized:
956     fprintf(stderr, "%s: unrecognized sender sampling: h=%s\n",
957     progname, curparams.hemis);
958     return(-1);
959     }
960    
961     /* Compute normal and area for polygon */
962     static int
963     finish_polygon(SURF *p)
964     {
965     const int nv = p->nfargs / 3;
966     FVECT e1, e2, vc;
967     int i;
968    
969     memset(p->snrm, 0, sizeof(FVECT));
970     VSUB(e1, p->farg+3, p->farg);
971     for (i = 2; i < nv; i++) {
972     VSUB(e2, p->farg+3*i, p->farg);
973     VCROSS(vc, e1, e2);
974     p->snrm[0] += vc[0];
975     p->snrm[1] += vc[1];
976     p->snrm[2] += vc[2];
977     VCOPY(e1, e2);
978     }
979     p->area = normalize(p->snrm)*0.5;
980     return(p->area > FTINY);
981     }
982    
983     /* Add a surface to our current parameters */
984     static void
985     add_surface(int st, const char *oname, FILE *fp)
986     {
987     SURF *snew;
988     int n;
989     /* get floating-point arguments */
990     if (!fscanf(fp, "%d", &n)) return;
991     while (n-- > 0) fscanf(fp, "%*s");
992     if (!fscanf(fp, "%d", &n)) return;
993     while (n-- > 0) fscanf(fp, "%*d");
994     if (!fscanf(fp, "%d", &n) || n <= 0) return;
995     snew = (SURF *)malloc(sizeof(SURF) + sizeof(double)*(n-1));
996     if (snew == NULL) {
997     fputs(progname, stderr);
998     fputs(": out of memory!\n", stderr);
999     exit(1);
1000     }
1001     strncpy(snew->sname, oname, sizeof(snew->sname)-1);
1002     snew->sname[sizeof(snew->sname)-1] = '\0';
1003     snew->styp = st;
1004     snew->priv = NULL;
1005     snew->nfargs = n;
1006     for (n = 0; n < snew->nfargs; n++)
1007     if (fscanf(fp, "%lf", &snew->farg[n]) != 1) {
1008     fprintf(stderr, "%s: error reading arguments for '%s'\n",
1009     progname, oname);
1010     exit(1);
1011     }
1012     switch (st) {
1013     case ST_RING:
1014     if (snew->nfargs != 8)
1015     goto badcount;
1016     VCOPY(snew->snrm, snew->farg+3);
1017     if (normalize(snew->snrm) == 0)
1018     goto badnorm;
1019     if (snew->farg[7] < snew->farg[6]) {
1020     double t = snew->farg[7];
1021     snew->farg[7] = snew->farg[6];
1022     snew->farg[6] = t;
1023     }
1024     snew->area = PI*(snew->farg[7]*snew->farg[7] -
1025     snew->farg[6]*snew->farg[6]);
1026     break;
1027     case ST_POLY:
1028     if (snew->nfargs < 9 || snew->nfargs % 3)
1029     goto badcount;
1030     finish_polygon(snew);
1031     break;
1032     case ST_SOURCE:
1033     if (snew->nfargs != 4)
1034     goto badcount;
1035 greg 2.8 for (n = 3; n--; ) /* need to reverse "normal" */
1036     snew->snrm[n] = -snew->farg[n];
1037 greg 2.1 if (normalize(snew->snrm) == 0)
1038     goto badnorm;
1039 greg 2.4 snew->area = sin((PI/180./2.)*snew->farg[3]);
1040     snew->area *= PI*snew->area;
1041 greg 2.1 break;
1042     }
1043 greg 2.16 if ((snew->area <= FTINY) & (verbose >= 0)) {
1044 greg 2.1 fprintf(stderr, "%s: warning - zero area for surface '%s'\n",
1045     progname, oname);
1046     free(snew);
1047     return;
1048     }
1049     VSUM(curparams.nrm, curparams.nrm, snew->snrm, snew->area);
1050     snew->next = curparams.slist;
1051     curparams.slist = snew;
1052     curparams.nsurfs++;
1053     return;
1054     badcount:
1055 greg 2.7 fprintf(stderr, "%s: bad argument count for surface element '%s'\n",
1056 greg 2.1 progname, oname);
1057     exit(1);
1058     badnorm:
1059 greg 2.7 fprintf(stderr, "%s: bad orientation for surface element '%s'\n",
1060 greg 2.1 progname, oname);
1061     exit(1);
1062     }
1063    
1064     /* Parse a receiver object (look for modifiers to add to rcontrib command) */
1065     static int
1066     add_recv_object(FILE *fp)
1067     {
1068     int st;
1069     char thismod[128], otype[32], oname[128];
1070     int n;
1071    
1072     if (fscanf(fp, "%s %s %s", thismod, otype, oname) != 3)
1073     return(0); /* must have hit EOF! */
1074     if (!strcmp(otype, "alias")) {
1075     fscanf(fp, "%*s"); /* skip alias */
1076     return(0);
1077     }
1078     /* is it a new receiver? */
1079     if ((st = surf_type(otype)) != ST_NONE) {
1080     if (curparams.slist != NULL && (st == ST_SOURCE) ^
1081     (curparams.slist->styp == ST_SOURCE)) {
1082     fputs(progname, stderr);
1083     fputs(": cannot mix source/non-source receivers!\n", stderr);
1084     return(-1);
1085     }
1086     if (strcmp(thismod, curmod)) {
1087     if (curmod[0]) { /* output last receiver? */
1088     finish_receiver();
1089     clear_params(&curparams, 1);
1090     }
1091 greg 2.6 parse_params(&curparams, newparams);
1092     newparams[0] = '\0';
1093 greg 2.1 strcpy(curmod, thismod);
1094     }
1095     add_surface(st, oname, fp); /* read & store surface */
1096     return(1);
1097     }
1098     /* else skip arguments */
1099 greg 2.8 if (!fscanf(fp, "%d", &n)) return(0);
1100 greg 2.1 while (n-- > 0) fscanf(fp, "%*s");
1101 greg 2.10 if (!fscanf(fp, "%d", &n)) return(0);
1102 greg 2.1 while (n-- > 0) fscanf(fp, "%*d");
1103 greg 2.10 if (!fscanf(fp, "%d", &n)) return(0);
1104 greg 2.1 while (n-- > 0) fscanf(fp, "%*f");
1105     return(0);
1106     }
1107    
1108     /* Parse a sender object */
1109     static int
1110     add_send_object(FILE *fp)
1111     {
1112     int st;
1113 greg 2.24 char thismod[128], otype[32], oname[128];
1114 greg 2.1 int n;
1115    
1116 greg 2.25 if (fscanf(fp, "%s %s %s", thismod, otype, oname) != 3)
1117 greg 2.1 return(0); /* must have hit EOF! */
1118     if (!strcmp(otype, "alias")) {
1119     fscanf(fp, "%*s"); /* skip alias */
1120     return(0);
1121     }
1122     /* is it a new surface? */
1123     if ((st = surf_type(otype)) != ST_NONE) {
1124     if (st == ST_SOURCE) {
1125     fputs(progname, stderr);
1126     fputs(": cannot use source as a sender!\n", stderr);
1127     return(-1);
1128     }
1129 greg 2.24 if (strcmp(thismod, curmod)) {
1130     if (curmod[0]) {
1131     fputs(progname, stderr);
1132     fputs(": warning - multiple modifiers in sender\n",
1133     stderr);
1134     }
1135     strcpy(curmod, thismod);
1136     }
1137 greg 2.6 parse_params(&curparams, newparams);
1138     newparams[0] = '\0';
1139 greg 2.1 add_surface(st, oname, fp); /* read & store surface */
1140     return(0);
1141     }
1142     /* else skip arguments */
1143 greg 2.8 if (!fscanf(fp, "%d", &n)) return(0);
1144 greg 2.1 while (n-- > 0) fscanf(fp, "%*s");
1145 greg 2.10 if (!fscanf(fp, "%d", &n)) return(0);
1146 greg 2.1 while (n-- > 0) fscanf(fp, "%*d");
1147 greg 2.10 if (!fscanf(fp, "%d", &n)) return(0);
1148 greg 2.1 while (n-- > 0) fscanf(fp, "%*f");
1149     return(0);
1150     }
1151    
1152     /* Load a Radiance scene using the given callback function for objects */
1153     static int
1154     load_scene(const char *inspec, int (*ocb)(FILE *))
1155     {
1156     int rv = 0;
1157     char inpbuf[1024];
1158     FILE *fp;
1159     int c;
1160    
1161     if (*inspec == '!')
1162     fp = popen(inspec+1, "r");
1163     else
1164     fp = fopen(inspec, "r");
1165     if (fp == NULL) {
1166     fprintf(stderr, "%s: cannot load '%s'\n", progname, inspec);
1167     return(-1);
1168     }
1169     while ((c = getc(fp)) != EOF) { /* load receiver data */
1170     if (isspace(c)) /* skip leading white space */
1171     continue;
1172     if (c == '!') { /* read from a new command */
1173     inpbuf[0] = c;
1174     if (fgetline(inpbuf+1, sizeof(inpbuf)-1, fp) != NULL) {
1175     if ((c = load_scene(inpbuf, ocb)) < 0)
1176     return(c);
1177     rv += c;
1178     }
1179     continue;
1180     }
1181     if (c == '#') { /* parameters/comment */
1182 greg 2.5 if ((c = getc(fp)) == EOF || ungetc(c,fp) == EOF)
1183     break;
1184     if (!isspace(c) && fscanf(fp, "%s", inpbuf) == 1 &&
1185 greg 2.1 !strcmp(inpbuf, PARAMSTART)) {
1186     if (fgets(inpbuf, sizeof(inpbuf), fp) != NULL)
1187 greg 2.6 strcat(newparams, inpbuf);
1188 greg 2.1 continue;
1189     }
1190 greg 2.10 while ((c = getc(fp)) != EOF && c != '\n')
1191 greg 2.1 ; /* else skipping comment */
1192     continue;
1193     }
1194     ungetc(c, fp); /* else check object for receiver */
1195     c = (*ocb)(fp);
1196     if (c < 0)
1197     return(c);
1198     rv += c;
1199     }
1200     /* close our input stream */
1201     c = (*inspec == '!') ? pclose(fp) : fclose(fp);
1202     if (c != 0) {
1203     fprintf(stderr, "%s: error loading '%s'\n", progname, inspec);
1204     return(-1);
1205     }
1206     return(rv);
1207     }
1208    
1209     /* Get command arguments and run program */
1210     int
1211     main(int argc, char *argv[])
1212     {
1213     char fmtopt[6] = "-faa"; /* default output is ASCII */
1214 greg 2.4 char *xrs=NULL, *yrs=NULL, *ldopt=NULL;
1215 greg 2.11 char *iropt = NULL;
1216 greg 2.1 char *sendfn;
1217     char sampcntbuf[32], nsbinbuf[32];
1218     FILE *rcfp;
1219     int nsbins;
1220     int a, i;
1221     /* screen rcontrib options */
1222     progname = argv[0];
1223 greg 2.14 for (a = 1; a < argc-2; a++) {
1224     int na;
1225     /* check for argument expansion */
1226     while ((na = expandarg(&argc, &argv, a)) > 0)
1227     ;
1228     if (na < 0) {
1229     fprintf(stderr, "%s: cannot expand '%s'\n",
1230     progname, argv[a]);
1231     return(1);
1232     }
1233     if (argv[a][0] != '-' || !argv[a][1])
1234     break;
1235     na = 1;
1236     switch (argv[a][1]) { /* !! Keep consistent !! */
1237 greg 2.1 case 'v': /* verbose mode */
1238 greg 2.16 verbose = 1;
1239 greg 2.1 na = 0;
1240     continue;
1241     case 'f': /* special case for -fo, -ff, etc. */
1242     switch (argv[a][2]) {
1243     case '\0': /* cal file */
1244     goto userr;
1245     case 'o': /* force output */
1246     goto userr;
1247     case 'a': /* output format */
1248     case 'f':
1249     case 'd':
1250     case 'c':
1251 greg 2.6 if (!(fmtopt[3] = argv[a][3]))
1252     fmtopt[3] = argv[a][2];
1253     fmtopt[2] = argv[a][2];
1254 greg 2.1 na = 0;
1255     continue; /* will pass later */
1256     default:
1257     goto userr;
1258     }
1259     break;
1260 greg 2.4 case 'x': /* x-resolution */
1261     xrs = argv[++a];
1262     na = 0;
1263     continue;
1264     case 'y': /* y-resolution */
1265     yrs = argv[++a];
1266     na = 0;
1267     continue;
1268 greg 2.1 case 'c': /* number of samples */
1269 greg 2.6 sampcnt = atoi(argv[++a]);
1270 greg 2.1 if (sampcnt <= 0)
1271     goto userr;
1272     na = 0; /* we re-add this later */
1273     continue;
1274 greg 2.8 case 'I': /* only for pass-through mode */
1275 greg 2.11 case 'i':
1276     iropt = argv[a];
1277 greg 2.8 na = 0;
1278     continue;
1279 greg 2.16 case 'w': /* options without arguments */
1280     if (argv[a][2] != '+') verbose = -1;
1281     case 'V':
1282 greg 2.1 case 'u':
1283     case 'h':
1284 greg 2.4 case 'r':
1285 greg 2.1 break;
1286     case 'n': /* options with 1 argument */
1287     case 's':
1288     case 'o':
1289     na = 2;
1290     break;
1291     case 'b': /* special case */
1292     if (argv[a][2] != 'v') goto userr;
1293     break;
1294     case 'l': /* special case */
1295 greg 2.4 if (argv[a][2] == 'd') {
1296     ldopt = argv[a];
1297     na = 0;
1298     continue;
1299     }
1300 greg 2.1 na = 2;
1301     break;
1302     case 'd': /* special case */
1303     if (argv[a][2] != 'v') na = 2;
1304     break;
1305     case 'a': /* special case */
1306 greg 2.28 if (argv[a][2] == 'p') {
1307     na = 2; /* photon map [+ bandwidth(s)] */
1308     if (a < argc-3 && atoi(argv[a+1]) > 0)
1309     na += 1 + (a < argc-4 && atoi(argv[a+2]) > 0);
1310     } else
1311     na = (argv[a][2] == 'v') ? 4 : 2;
1312 greg 2.1 break;
1313     case 'm': /* special case */
1314     if (!argv[a][2]) goto userr;
1315     na = (argv[a][2] == 'e') | (argv[a][2] == 'a') ? 4 : 2;
1316     break;
1317     default: /* anything else is verbotten */
1318     goto userr;
1319     }
1320     if (na <= 0) continue;
1321     CHECKARGC(na); /* pass on option */
1322     rcarg[nrcargs++] = argv[a];
1323     while (--na) /* + arguments if any */
1324     rcarg[nrcargs++] = argv[++a];
1325     }
1326     if (a > argc-2)
1327     goto userr; /* check at end of options */
1328     sendfn = argv[a++]; /* assign sender & receiver inputs */
1329     if (sendfn[0] == '-') { /* user wants pass-through mode? */
1330     if (sendfn[1]) goto userr;
1331     sendfn = NULL;
1332 greg 2.11 if (iropt) {
1333 greg 2.8 CHECKARGC(1);
1334 greg 2.11 rcarg[nrcargs++] = iropt;
1335 greg 2.8 }
1336 greg 2.4 if (xrs) {
1337     CHECKARGC(2);
1338     rcarg[nrcargs++] = "-x";
1339     rcarg[nrcargs++] = xrs;
1340     }
1341     if (yrs) {
1342     CHECKARGC(2);
1343     rcarg[nrcargs++] = "-y";
1344     rcarg[nrcargs++] = yrs;
1345     }
1346     if (ldopt) {
1347     CHECKARGC(1);
1348     rcarg[nrcargs++] = ldopt;
1349     }
1350 greg 2.1 if (sampcnt <= 0) sampcnt = 1;
1351 greg 2.8 } else { /* else in sampling mode */
1352 greg 2.11 if (iropt) {
1353 greg 2.8 fputs(progname, stderr);
1354 greg 2.11 fputs(": -i, -I supported for pass-through only\n", stderr);
1355 greg 2.8 return(1);
1356     }
1357 greg 2.6 fmtopt[2] = (sizeof(RREAL)==sizeof(double)) ? 'd' : 'f';
1358 greg 2.1 if (sampcnt <= 0) sampcnt = 10000;
1359     }
1360     sprintf(sampcntbuf, "%d", sampcnt);
1361     CHECKARGC(3); /* add our format & sample count */
1362     rcarg[nrcargs++] = fmtopt;
1363     rcarg[nrcargs++] = "-c";
1364     rcarg[nrcargs++] = sampcntbuf;
1365     /* add receiver arguments to rcontrib */
1366     if (load_scene(argv[a], add_recv_object) < 0)
1367     return(1);
1368     finish_receiver();
1369     if (sendfn == NULL) { /* pass-through mode? */
1370     CHECKARGC(1); /* add octree */
1371     rcarg[nrcargs++] = oconv_command(argc-a, argv+a);
1372     rcarg[nrcargs] = NULL;
1373     return(my_exec(rcarg)); /* rcontrib does everything */
1374     }
1375     clear_params(&curparams, 0); /* else load sender surface & params */
1376 greg 2.24 curmod[0] = '\0';
1377 greg 2.1 if (load_scene(sendfn, add_send_object) < 0)
1378     return(1);
1379     if ((nsbins = prepare_sampler()) <= 0)
1380     return(1);
1381     CHECKARGC(3); /* add row count and octree */
1382     rcarg[nrcargs++] = "-y";
1383     sprintf(nsbinbuf, "%d", nsbins);
1384     rcarg[nrcargs++] = nsbinbuf;
1385     rcarg[nrcargs++] = oconv_command(argc-a, argv+a);
1386     rcarg[nrcargs] = NULL;
1387     /* open pipe to rcontrib process */
1388     if ((rcfp = popen_arglist(rcarg, "w")) == NULL)
1389     return(1);
1390     SET_FILE_BINARY(rcfp);
1391     #ifdef getc_unlocked
1392     flockfile(rcfp);
1393     #endif
1394 greg 2.22 if (verbose > 0) {
1395 greg 2.1 fprintf(stderr, "%s: sampling %d directions", progname, nsbins);
1396     if (curparams.nsurfs > 1)
1397 greg 2.7 fprintf(stderr, " (%d elements)\n", curparams.nsurfs);
1398 greg 2.1 else
1399     fputc('\n', stderr);
1400     }
1401     for (i = 0; i < nsbins; i++) /* send rcontrib ray samples */
1402     if (!(*curparams.sample_basis)(&curparams, i, rcfp))
1403     return(1);
1404 greg 2.23 return(pclose(rcfp) < 0); /* all finished! */
1405 greg 2.1 userr:
1406     if (a < argc-2)
1407     fprintf(stderr, "%s: unsupported option '%s'", progname, argv[a]);
1408 greg 2.13 fprintf(stderr, "Usage: %s [-v][rcontrib options] sender.rad receiver.rad [-i system.oct] [system.rad ..]\n",
1409 greg 2.1 progname);
1410     return(1);
1411     }