ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rfluxmtx.c
Revision: 2.5
Committed: Tue Jul 22 23:21:56 2014 UTC (9 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +7 -4 lines
Log Message:
Fixed some bugs... more to go

File Contents

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