ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.7
Committed: Sat Dec 12 19:01:00 2009 UTC (14 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0
Changes since 2.6: +12 -1 lines
Log Message:
Added -n option to rtrace and moved quit() funciton out of raypcalls

File Contents

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