ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.5
Committed: Fri Apr 11 22:06:04 2008 UTC (15 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R9
Changes since 2.4: +8 -4 lines
Log Message:
Changed default view and added handling of fore clipping distance

File Contents

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