ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.27
Committed: Wed Apr 1 01:38:36 2015 UTC (9 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.26: +2 -2 lines
Log Message:
Fixes for Klems & tensor tree coordinate orientations

File Contents

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