ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.3
Committed: Fri Feb 22 18:04:02 2008 UTC (16 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +23 -19 lines
Log Message:
Fixed options processing

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.3 static const char RCSid[] = "$Id: rsensor.c,v 2.2 2008/02/21 20:18:25 greg Exp $";
3 greg 2.1 #endif
4    
5     /*
6     * Compute sensor signal based on spatial sensitivity.
7     *
8     * Created Feb 2008 for Architectural Energy Corp.
9     */
10    
11     #include "ray.h"
12     #include "source.h"
13     #include "view.h"
14     #include "random.h"
15    
16     #define DEGREE (PI/180.)
17    
18     #define MAXNT 180 /* maximum number of theta divisions */
19     #define MAXNP 360 /* maximum number of phi divisions */
20    
21     extern char *progname; /* global argv[0] */
22     extern int nowarn; /* don't report warnings? */
23    
24     /* current sensor's perspective */
25     VIEW ourview = STDVIEW;
26    
27     unsigned long nsamps = 10000; /* desired number of initial samples */
28     unsigned long nssamps = 9000; /* number of super-samples */
29 greg 2.2 int ndsamps = 32; /* number of direct samples */
30 greg 2.1 int nprocs = 1; /* number of rendering processes */
31    
32     float *sensor = NULL; /* current sensor data */
33     int sntp[2]; /* number of sensor theta and phi angles */
34     float maxtheta; /* maximum theta value for this sensor */
35     float tvals[MAXNT+1]; /* theta values (1-D table of 1-cos(t)) */
36     float *pvals = NULL; /* phi values (2-D table in radians) */
37     int ntheta = 0; /* polar angle divisions */
38     int nphi = 0; /* azimuthal angle divisions */
39     double gscale = 1.; /* global scaling value */
40    
41 greg 2.2 #define s_theta(t) sensor[(t+1)*(sntp[1]+1)]
42     #define s_phi(p) sensor[(p)+1]
43     #define s_val(t,p) sensor[(p)+1+(t+1)*(sntp[1]+1)]
44    
45 greg 2.1 static void comp_sensor(char *sfile);
46    
47     static void
48 greg 2.3 over_options() /* overriding options */
49 greg 2.1 {
50 greg 2.3 directvis = (ndsamps <= 0);
51     do_irrad = 0;
52     }
53    
54     static void
55     print_defaults() /* print out default parameters */
56     {
57     over_options();
58 greg 2.1 printf("-n %-9d\t\t\t# number of processes\n", nprocs);
59     printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
60     /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
61     printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
62     printf("-vp %f %f %f\t# view point\n",
63     ourview.vp[0], ourview.vp[1], ourview.vp[2]);
64     printf("-vd %f %f %f\t# view direction\n",
65     ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]);
66     printf("-vu %f %f %f\t# view up\n",
67     ourview.vup[0], ourview.vup[1], ourview.vup[2]);
68     printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore);
69     print_rdefaults();
70     }
71    
72     int
73     main(
74     int argc,
75     char *argv[]
76     )
77     {
78     int doheader = 1;
79 greg 2.3 int optwarn = 0;
80 greg 2.1 int i, rval;
81    
82     progname = argv[0];
83     /* set up rendering defaults */
84 greg 2.2 rand_samp = 1;
85     dstrsrc = 0.5;
86     srcsizerat = 0.1;
87 greg 2.1 directrelay = 3;
88     ambounce = 1;
89 greg 2.2 maxdepth = -10;
90 greg 2.1 /* get options from command line */
91 greg 2.3 for (i = 1; i < argc; i++) {
92 greg 2.1 while ((rval = expandarg(&argc, &argv, i)) > 0)
93     ;
94     if (rval < 0) {
95     sprintf(errmsg, "cannot expand '%s'", argv[i]);
96     error(SYSTEM, errmsg);
97     }
98 greg 2.3 if (argv[i][0] != '-') {
99     if (i >= argc-1)
100     break; /* final octree argument */
101 greg 2.1 if (!ray_pnprocs) {
102 greg 2.3 over_options();
103 greg 2.1 if (doheader) { /* print header */
104     printargs(argc, argv, stdout);
105     fputformat("ascii", stdout);
106     putchar('\n');
107     }
108     /* start process(es) */
109     ray_pinit(argv[argc-1], nprocs);
110     }
111 greg 2.3 comp_sensor(argv[i]); /* process a sensor file */
112 greg 2.1 continue;
113     }
114     if (argv[i][1] == 'r') { /* sampling options */
115     if (argv[i][2] == 'd')
116     nsamps = atol(argv[++i]);
117     else if (argv[i][2] == 's')
118     nssamps = atol(argv[++i]);
119     else {
120     sprintf(errmsg, "bad option at '%s'", argv[i]);
121     error(USER, errmsg);
122     }
123     continue;
124     }
125     /* direct component samples */
126     if (argv[i][1] == 'd' && argv[i][2] == 'n') {
127     ndsamps = atoi(argv[++i]);
128     continue;
129     }
130     if (argv[i][1] == 'v') { /* next sensor view */
131     if (argv[i][2] == 'f') {
132     rval = viewfile(argv[++i], &ourview, NULL);
133     if (rval < 0) {
134     sprintf(errmsg,
135     "cannot open view file \"%s\"",
136     argv[i]);
137     error(SYSTEM, errmsg);
138     } else if (rval == 0) {
139     sprintf(errmsg,
140     "bad view file \"%s\"",
141     argv[i]);
142     error(USER, errmsg);
143     }
144     continue;
145     }
146     rval = getviewopt(&ourview, argc-i, argv+i);
147     if (rval >= 0) {
148     i += rval;
149     continue;
150     }
151     sprintf(errmsg, "bad view option at '%s'", argv[i]);
152     error(USER, errmsg);
153     }
154 greg 2.3 if (!strcmp(argv[i], "-w")) { /* toggle warnings */
155     nowarn = !nowarn;
156 greg 2.1 continue;
157     }
158     if (ray_pnprocs) {
159 greg 2.3 if (!optwarn++)
160     error(WARNING,
161 greg 2.1 "rendering options should appear before first sensor");
162     } else if (!strcmp(argv[i], "-defaults")) {
163     print_defaults();
164     return(0);
165     }
166     if (argv[i][1] == 'h') { /* header toggle */
167     doheader = !doheader;
168     continue;
169     }
170     if (!strcmp(argv[i], "-n")) { /* number of processes */
171     nprocs = atoi(argv[++i]);
172     if (nprocs <= 0)
173     error(USER, "illegal number of processes");
174     continue;
175     }
176     rval = getrenderopt(argc-i, argv+i);
177     if (rval < 0) {
178     sprintf(errmsg, "bad render option at '%s'", argv[i]);
179     error(USER, errmsg);
180     }
181     i += rval;
182     }
183 greg 2.3 if (!ray_pnprocs)
184     error(USER, i<argc ? "missing sensor file" : "missing octree");
185 greg 2.1 quit(0);
186     }
187    
188     /* Load sensor sensitivities (first row and column are angles) */
189     static float *
190     load_sensor(
191     int ntp[2],
192     char *sfile
193     )
194     {
195     char linebuf[8192];
196     int nelem = 1000;
197     float *sarr = (float *)malloc(sizeof(float)*nelem);
198     FILE *fp;
199     char *cp;
200     int i;
201    
202     fp = frlibopen(sfile);
203     if (fp == NULL) {
204     sprintf(errmsg, "cannot open sensor file '%s'", sfile);
205     error(SYSTEM, errmsg);
206     }
207     fgets(linebuf, sizeof(linebuf), fp);
208     if (!strncmp(linebuf, "Elevation ", 10))
209     fgets(linebuf, sizeof(linebuf), fp);
210     /* get phi values */
211     sarr[0] = .0f;
212     if (strncmp(linebuf, "degrees", 7)) {
213     sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
214     error(USER, errmsg);
215     }
216     cp = sskip(linebuf);
217     ntp[1] = 0;
218     for ( ; ; ) {
219     sarr[ntp[1]+1] = atof(cp);
220     cp = fskip(cp);
221     if (cp == NULL)
222     break;
223     ++ntp[1];
224     }
225     ntp[0] = 0; /* get thetas + data */
226     while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
227     ++ntp[0];
228     if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
229     nelem += (nelem>>2) + ntp[1];
230     sarr = (float *)realloc((void *)sarr,
231     sizeof(float)*nelem);
232     if (sarr == NULL)
233     error(SYSTEM, "out of memory in load_sensor()");
234     }
235     cp = linebuf;
236     i = ntp[0]*(ntp[1]+1);
237     for ( ; ; ) {
238     sarr[i] = atof(cp);
239     cp = fskip(cp);
240     if (cp == NULL)
241     break;
242     ++i;
243     }
244     if (i == ntp[0]*(ntp[1]+1))
245     break;
246     if (i != (ntp[0]+1)*(ntp[1]+1)) {
247     sprintf(errmsg,
248     "bad column count near line %d in sensor file '%s'",
249     ntp[0]+1, sfile);
250     error(USER, errmsg);
251     }
252     }
253     nelem = i;
254     fclose(fp);
255     errmsg[0] = '\0'; /* sanity checks */
256     if (ntp[0] <= 0)
257     sprintf(errmsg, "no data in sensor file '%s'", sfile);
258     else if (fabs(sarr[ntp[1]+1]) > FTINY)
259     sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
260     sfile);
261     else if (fabs(sarr[1]) > FTINY)
262     sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
263     sfile);
264     else if (sarr[ntp[1]] <= FTINY)
265     sprintf(errmsg,
266     "maximum phi must be positive in sensor file '%s'",
267     sfile);
268     else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
269     sprintf(errmsg,
270     "maximum theta must be positive in sensor file '%s'",
271     sfile);
272     if (errmsg[0])
273     error(USER, errmsg);
274     return((float *)realloc((void *)sarr, sizeof(float)*nelem));
275     }
276    
277     /* Initialize probability table */
278     static void
279     init_ptable(
280     char *sfile
281     )
282     {
283     int samptot = nsamps;
284     float *rowp, *rowp1;
285     double rowsum[MAXNT], rowomega[MAXNT];
286     double thdiv[MAXNT+1], phdiv[MAXNP+1];
287     double tsize, psize;
288     double prob, frac, frac1;
289     int i, j, t, p;
290     /* free old table */
291     if (sensor != NULL)
292     free((void *)sensor);
293     if (pvals != NULL)
294     free((void *)pvals);
295     if (sfile == NULL || !*sfile) {
296 greg 2.2 sensor = NULL;
297     sntp[0] = sntp[1] = 0;
298 greg 2.1 pvals = NULL;
299     ntheta = nphi = 0;
300     return;
301     }
302     /* load sensor table */
303     sensor = load_sensor(sntp, sfile);
304     if (sntp[0] > MAXNT) {
305     sprintf(errmsg, "Too many theta rows in sensor file '%s'",
306     sfile);
307     error(INTERNAL, errmsg);
308     }
309     if (sntp[1] > MAXNP) {
310     sprintf(errmsg, "Too many phi columns in sensor file '%s'",
311     sfile);
312     error(INTERNAL, errmsg);
313     }
314     /* compute boundary angles */
315 greg 2.2 maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
316 greg 2.1 thdiv[0] = .0;
317     for (t = 1; t < sntp[0]; t++)
318 greg 2.2 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
319 greg 2.1 thdiv[sntp[0]] = maxtheta*DEGREE;
320     phdiv[0] = .0;
321     for (p = 1; p < sntp[1]; p++)
322 greg 2.2 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
323 greg 2.1 phdiv[sntp[1]] = 2.*PI;
324     /* size our table */
325     tsize = 1. - cos(maxtheta*DEGREE);
326 greg 2.2 psize = PI*tsize/(maxtheta*DEGREE);
327 greg 2.1 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
328     samptot = sntp[0]*sntp[1];
329 greg 2.2 ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
330 greg 2.1 if (ntheta > MAXNT)
331     ntheta = MAXNT;
332     nphi = samptot/ntheta;
333     pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
334     if (pvals == NULL)
335     error(SYSTEM, "out of memory in init_ptable()");
336     gscale = .0; /* compute our inverse table */
337     for (i = 0; i < sntp[0]; i++) {
338 greg 2.2 rowp = &s_val(i,0);
339 greg 2.1 rowsum[i] = 0.;
340     for (j = 0; j < sntp[1]; j++)
341     rowsum[i] += *rowp++;
342     rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
343     rowomega[i] *= 2.*PI / (double)sntp[1];
344     gscale += rowsum[i] * rowomega[i];
345     }
346 greg 2.2 for (i = 0; i < ntheta; i++) {
347 greg 2.1 prob = (double)i / (double)ntheta;
348     for (t = 0; t < sntp[0]; t++)
349     if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
350     break;
351     if (t >= sntp[0])
352     error(INTERNAL, "code error 1 in init_ptable()");
353     frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
354     tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
355     frac*cos(thdiv[t+1]) );
356 greg 2.2 /* offset b/c sensor values are centered */
357     if (t <= 0 || frac > 0.5)
358     frac -= 0.5;
359     else if (t >= sntp[0]-1 || frac < 0.5) {
360     frac += 0.5;
361     --t;
362     }
363 greg 2.1 pvals[i*(nphi+1)] = .0f;
364     for (j = 1; j < nphi; j++) {
365     prob = (double)j / (double)nphi;
366 greg 2.2 rowp = &s_val(t,0);
367     rowp1 = &s_val(t+1,0);
368 greg 2.1 for (p = 0; p < sntp[1]; p++) {
369 greg 2.2 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
370     frac*rowp1[p]/rowsum[t+1]) <= .0)
371 greg 2.1 break;
372     if (p >= sntp[1])
373     error(INTERNAL,
374     "code error 2 in init_ptable()");
375 greg 2.2 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
376     + frac*rowp1[p]/rowsum[t+1]);
377     if (p <= 0 || frac1 > 0.5)
378     frac1 -= 0.5;
379     else if (p >= sntp[1]-1 || frac1 < 0.5) {
380     frac1 += 0.5;
381     --p;
382     }
383 greg 2.1 pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
384     frac1*phdiv[p+1];
385     }
386     }
387     pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
388     }
389 greg 2.2 tvals[0] = .0f;
390 greg 2.1 tvals[ntheta] = (float)tsize;
391     }
392    
393     /* Get normalized direction from random variables in [0,1) range */
394     static void
395     get_direc(
396     FVECT dvec,
397     double x,
398     double y
399     )
400     {
401     double xfrac = x*ntheta;
402     int tndx = (int)xfrac;
403     double yfrac = y*nphi;
404     int pndx = (int)yfrac;
405     double rad, phi;
406     FVECT dv;
407     int i;
408    
409     xfrac -= (double)tndx;
410     yfrac -= (double)pndx;
411     pndx += tndx*(nphi+1);
412    
413     dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]);
414     rad = sqrt(1. - dv[2]*dv[2]);
415     phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1];
416     dv[0] = -rad*sin(phi);
417     dv[1] = rad*cos(phi);
418     for (i = 3; i--; )
419     dvec[i] = dv[0]*ourview.hvec[i] +
420     dv[1]*ourview.vvec[i] +
421     dv[2]*ourview.vdir[i] ;
422     }
423    
424     /* Get sensor value in the specified direction (normalized) */
425     static float
426     sens_val(
427     FVECT dvec
428     )
429     {
430     FVECT dv;
431     float theta, phi;
432     int t, p;
433    
434     dv[2] = DOT(dvec, ourview.vdir);
435     theta = (float)((1./DEGREE) * acos(dv[2]));
436     if (theta >= maxtheta)
437     return(.0f);
438     dv[0] = DOT(dvec, ourview.hvec);
439     dv[1] = DOT(dvec, ourview.vvec);
440     phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
441     while (phi < .0f) phi += 360.f;
442     t = (int)(theta/maxtheta * sntp[0]);
443     p = (int)(phi*(1./360.) * sntp[1]);
444     /* hack for non-uniform sensor grid */
445 greg 2.2 while (t+1 < sntp[0] && theta >= s_theta(t+1))
446 greg 2.1 ++t;
447 greg 2.2 while (t-1 >= 0 && theta <= s_theta(t-1))
448 greg 2.1 --t;
449 greg 2.2 while (p+1 < sntp[1] && phi >= s_phi(p+1))
450 greg 2.1 ++p;
451 greg 2.2 while (p-1 >= 0 && phi <= s_phi(p-1))
452 greg 2.1 --p;
453 greg 2.2 return(s_val(t,p));
454 greg 2.1 }
455    
456     /* Compute sensor output */
457     static void
458     comp_sensor(
459     char *sfile
460     )
461     {
462     int ndirs = dstrsrc > FTINY ? ndsamps :
463     ndsamps > 0 ? 1 : 0;
464     char *err;
465     int nt, np;
466     COLOR vsum;
467     RAY rr;
468     int i, j;
469     /* set view */
470     ourview.type = VT_ANG;
471     ourview.horiz = ourview.vert = 180.;
472     ourview.hoff = ourview.voff = .0;
473     err = setview(&ourview);
474     if (err != NULL)
475     error(USER, err);
476     /* assign probability table */
477     init_ptable(sfile);
478     /* do Monte Carlo sampling */
479     setcolor(vsum, .0f, .0f, .0f);
480     nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
481     np = nsamps/nt;
482     VCOPY(rr.rorg, ourview.vp);
483     rr.rmax = .0;
484     for (i = 0; i < nt; i++)
485 greg 2.2 for (j = 0; j < np; j++) {
486     get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
487 greg 2.1 rayorigin(&rr, PRIMARY, NULL, NULL);
488     if (ray_pqueue(&rr) == 1)
489     addcolor(vsum, rr.rcol);
490     }
491     /* finish MC calculation */
492     while (ray_presult(&rr, 0) > 0)
493     addcolor(vsum, rr.rcol);
494     scalecolor(vsum, gscale/(nt*np));
495     /* compute direct component */
496     for (i = ndirs; i-- > 0; ) {
497     SRCINDEX si;
498     initsrcindex(&si);
499     while (srcray(&rr, NULL, &si)) {
500     double d = sens_val(rr.rdir);
501     if (d <= FTINY)
502     continue;
503     d *= si.dom/ndirs;
504     scalecolor(rr.rcoef, d);
505     if (ray_pqueue(&rr) == 1) {
506     multcolor(rr.rcol, rr.rcoef);
507     addcolor(vsum, rr.rcol);
508     }
509     }
510     }
511     /* finish direct calculation */
512     while (ray_presult(&rr, 0) > 0) {
513     multcolor(rr.rcol, rr.rcoef);
514     addcolor(vsum, rr.rcol);
515     }
516     /* print our result */
517     printf("%.4e %.4e %.4e\n", colval(vsum,RED),
518     colval(vsum,GRN), colval(vsum,BLU));
519     }