ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.2
Committed: Thu Feb 21 20:18:25 2008 UTC (16 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +46 -27 lines
Log Message:
First tested version

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.2 static const char RCSid[] = "$Id: rsensor.c,v 2.1 2008/02/21 01:22:06 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     print_defaults()
49     {
50     printf("-n %-9d\t\t\t# number of processes\n", nprocs);
51     printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
52     /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
53     printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
54     printf("-vp %f %f %f\t# view point\n",
55     ourview.vp[0], ourview.vp[1], ourview.vp[2]);
56     printf("-vd %f %f %f\t# view direction\n",
57     ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]);
58     printf("-vu %f %f %f\t# view up\n",
59     ourview.vup[0], ourview.vup[1], ourview.vup[2]);
60     printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore);
61     print_rdefaults();
62     }
63    
64     int
65     main(
66     int argc,
67     char *argv[]
68     )
69     {
70     int doheader = 1;
71     int i, rval;
72    
73     progname = argv[0];
74     /* set up rendering defaults */
75 greg 2.2 rand_samp = 1;
76     dstrsrc = 0.5;
77     srcsizerat = 0.1;
78 greg 2.1 directrelay = 3;
79     ambounce = 1;
80 greg 2.2 maxdepth = -10;
81 greg 2.1 /* just asking defaults? */
82     if (argc == 2 && !strcmp(argv[1], "-defaults")) {
83     print_defaults();
84     return(0);
85     }
86     /* check octree */
87     if (argc < 2 || argv[argc-1][0] == '-')
88     error(USER, "missing octree argument");
89     /* get options from command line */
90     for (i = 1; i < argc-1; i++) {
91     while ((rval = expandarg(&argc, &argv, i)) > 0)
92     ;
93     if (rval < 0) {
94     sprintf(errmsg, "cannot expand '%s'", argv[i]);
95     error(SYSTEM, errmsg);
96     }
97     if (argv[i][0] != '-') { /* process a sensor file */
98     if (!ray_pnprocs) {
99     /* overriding options */
100     directvis = (ndsamps <= 0);
101     do_irrad = 0;
102     if (doheader) { /* print header */
103     printargs(argc, argv, stdout);
104     fputformat("ascii", stdout);
105     putchar('\n');
106     }
107     /* start process(es) */
108     ray_pinit(argv[argc-1], nprocs);
109     }
110     comp_sensor(argv[i]);
111     continue;
112     }
113     if (argv[i][1] == 'r') { /* sampling options */
114     if (argv[i][2] == 'd')
115     nsamps = atol(argv[++i]);
116     else if (argv[i][2] == 's')
117     nssamps = atol(argv[++i]);
118     else {
119     sprintf(errmsg, "bad option at '%s'", argv[i]);
120     error(USER, errmsg);
121     }
122     continue;
123     }
124     /* direct component samples */
125     if (argv[i][1] == 'd' && argv[i][2] == 'n') {
126     ndsamps = atoi(argv[++i]);
127     continue;
128     }
129     if (argv[i][1] == 'v') { /* next sensor view */
130     if (argv[i][2] == 'f') {
131     rval = viewfile(argv[++i], &ourview, NULL);
132     if (rval < 0) {
133     sprintf(errmsg,
134     "cannot open view file \"%s\"",
135     argv[i]);
136     error(SYSTEM, errmsg);
137     } else if (rval == 0) {
138     sprintf(errmsg,
139     "bad view file \"%s\"",
140     argv[i]);
141     error(USER, errmsg);
142     }
143     continue;
144     }
145     rval = getviewopt(&ourview, argc-i, argv+i);
146     if (rval >= 0) {
147     i += rval;
148     continue;
149     }
150     sprintf(errmsg, "bad view option at '%s'", argv[i]);
151     error(USER, errmsg);
152     }
153     if (!strcmp(argv[i], "-w")) { /* turn off warnings */
154     nowarn = 1;
155     continue;
156     }
157     if (ray_pnprocs) {
158     error(WARNING,
159     "rendering options should appear before first sensor");
160     } else if (!strcmp(argv[i], "-defaults")) {
161     print_defaults();
162     return(0);
163     }
164     if (argv[i][1] == 'h') { /* header toggle */
165     doheader = !doheader;
166     continue;
167     }
168     if (!strcmp(argv[i], "-n")) { /* number of processes */
169     nprocs = atoi(argv[++i]);
170     if (nprocs <= 0)
171     error(USER, "illegal number of processes");
172     continue;
173     }
174     rval = getrenderopt(argc-i, argv+i);
175     if (rval < 0) {
176     sprintf(errmsg, "bad render option at '%s'", argv[i]);
177     error(USER, errmsg);
178     }
179     i += rval;
180     }
181     quit(0);
182     }
183    
184     /* Load sensor sensitivities (first row and column are angles) */
185     static float *
186     load_sensor(
187     int ntp[2],
188     char *sfile
189     )
190     {
191     char linebuf[8192];
192     int nelem = 1000;
193     float *sarr = (float *)malloc(sizeof(float)*nelem);
194     FILE *fp;
195     char *cp;
196     int i;
197    
198     fp = frlibopen(sfile);
199     if (fp == NULL) {
200     sprintf(errmsg, "cannot open sensor file '%s'", sfile);
201     error(SYSTEM, errmsg);
202     }
203     fgets(linebuf, sizeof(linebuf), fp);
204     if (!strncmp(linebuf, "Elevation ", 10))
205     fgets(linebuf, sizeof(linebuf), fp);
206     /* get phi values */
207     sarr[0] = .0f;
208     if (strncmp(linebuf, "degrees", 7)) {
209     sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
210     error(USER, errmsg);
211     }
212     cp = sskip(linebuf);
213     ntp[1] = 0;
214     for ( ; ; ) {
215     sarr[ntp[1]+1] = atof(cp);
216     cp = fskip(cp);
217     if (cp == NULL)
218     break;
219     ++ntp[1];
220     }
221     ntp[0] = 0; /* get thetas + data */
222     while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
223     ++ntp[0];
224     if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
225     nelem += (nelem>>2) + ntp[1];
226     sarr = (float *)realloc((void *)sarr,
227     sizeof(float)*nelem);
228     if (sarr == NULL)
229     error(SYSTEM, "out of memory in load_sensor()");
230     }
231     cp = linebuf;
232     i = ntp[0]*(ntp[1]+1);
233     for ( ; ; ) {
234     sarr[i] = atof(cp);
235     cp = fskip(cp);
236     if (cp == NULL)
237     break;
238     ++i;
239     }
240     if (i == ntp[0]*(ntp[1]+1))
241     break;
242     if (i != (ntp[0]+1)*(ntp[1]+1)) {
243     sprintf(errmsg,
244     "bad column count near line %d in sensor file '%s'",
245     ntp[0]+1, sfile);
246     error(USER, errmsg);
247     }
248     }
249     nelem = i;
250     fclose(fp);
251     errmsg[0] = '\0'; /* sanity checks */
252     if (ntp[0] <= 0)
253     sprintf(errmsg, "no data in sensor file '%s'", sfile);
254     else if (fabs(sarr[ntp[1]+1]) > FTINY)
255     sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
256     sfile);
257     else if (fabs(sarr[1]) > FTINY)
258     sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
259     sfile);
260     else if (sarr[ntp[1]] <= FTINY)
261     sprintf(errmsg,
262     "maximum phi must be positive in sensor file '%s'",
263     sfile);
264     else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
265     sprintf(errmsg,
266     "maximum theta must be positive in sensor file '%s'",
267     sfile);
268     if (errmsg[0])
269     error(USER, errmsg);
270     return((float *)realloc((void *)sarr, sizeof(float)*nelem));
271     }
272    
273     /* Initialize probability table */
274     static void
275     init_ptable(
276     char *sfile
277     )
278     {
279     int samptot = nsamps;
280     float *rowp, *rowp1;
281     double rowsum[MAXNT], rowomega[MAXNT];
282     double thdiv[MAXNT+1], phdiv[MAXNP+1];
283     double tsize, psize;
284     double prob, frac, frac1;
285     int i, j, t, p;
286     /* free old table */
287     if (sensor != NULL)
288     free((void *)sensor);
289     if (pvals != NULL)
290     free((void *)pvals);
291     if (sfile == NULL || !*sfile) {
292 greg 2.2 sensor = NULL;
293     sntp[0] = sntp[1] = 0;
294 greg 2.1 pvals = NULL;
295     ntheta = nphi = 0;
296     return;
297     }
298     /* load sensor table */
299     sensor = load_sensor(sntp, sfile);
300     if (sntp[0] > MAXNT) {
301     sprintf(errmsg, "Too many theta rows in sensor file '%s'",
302     sfile);
303     error(INTERNAL, errmsg);
304     }
305     if (sntp[1] > MAXNP) {
306     sprintf(errmsg, "Too many phi columns in sensor file '%s'",
307     sfile);
308     error(INTERNAL, errmsg);
309     }
310     /* compute boundary angles */
311 greg 2.2 maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
312 greg 2.1 thdiv[0] = .0;
313     for (t = 1; t < sntp[0]; t++)
314 greg 2.2 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
315 greg 2.1 thdiv[sntp[0]] = maxtheta*DEGREE;
316     phdiv[0] = .0;
317     for (p = 1; p < sntp[1]; p++)
318 greg 2.2 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
319 greg 2.1 phdiv[sntp[1]] = 2.*PI;
320     /* size our table */
321     tsize = 1. - cos(maxtheta*DEGREE);
322 greg 2.2 psize = PI*tsize/(maxtheta*DEGREE);
323 greg 2.1 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
324     samptot = sntp[0]*sntp[1];
325 greg 2.2 ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
326 greg 2.1 if (ntheta > MAXNT)
327     ntheta = MAXNT;
328     nphi = samptot/ntheta;
329     pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
330     if (pvals == NULL)
331     error(SYSTEM, "out of memory in init_ptable()");
332     gscale = .0; /* compute our inverse table */
333     for (i = 0; i < sntp[0]; i++) {
334 greg 2.2 rowp = &s_val(i,0);
335 greg 2.1 rowsum[i] = 0.;
336     for (j = 0; j < sntp[1]; j++)
337     rowsum[i] += *rowp++;
338     rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
339     rowomega[i] *= 2.*PI / (double)sntp[1];
340     gscale += rowsum[i] * rowomega[i];
341     }
342 greg 2.2 for (i = 0; i < ntheta; i++) {
343 greg 2.1 prob = (double)i / (double)ntheta;
344     for (t = 0; t < sntp[0]; t++)
345     if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
346     break;
347     if (t >= sntp[0])
348     error(INTERNAL, "code error 1 in init_ptable()");
349     frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
350     tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
351     frac*cos(thdiv[t+1]) );
352 greg 2.2 /* offset b/c sensor values are centered */
353     if (t <= 0 || frac > 0.5)
354     frac -= 0.5;
355     else if (t >= sntp[0]-1 || frac < 0.5) {
356     frac += 0.5;
357     --t;
358     }
359 greg 2.1 pvals[i*(nphi+1)] = .0f;
360     for (j = 1; j < nphi; j++) {
361     prob = (double)j / (double)nphi;
362 greg 2.2 rowp = &s_val(t,0);
363     rowp1 = &s_val(t+1,0);
364 greg 2.1 for (p = 0; p < sntp[1]; p++) {
365 greg 2.2 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
366     frac*rowp1[p]/rowsum[t+1]) <= .0)
367 greg 2.1 break;
368     if (p >= sntp[1])
369     error(INTERNAL,
370     "code error 2 in init_ptable()");
371 greg 2.2 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
372     + frac*rowp1[p]/rowsum[t+1]);
373     if (p <= 0 || frac1 > 0.5)
374     frac1 -= 0.5;
375     else if (p >= sntp[1]-1 || frac1 < 0.5) {
376     frac1 += 0.5;
377     --p;
378     }
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     VCOPY(rr.rorg, ourview.vp);
479     rr.rmax = .0;
480     for (i = 0; i < nt; i++)
481 greg 2.2 for (j = 0; j < np; j++) {
482     get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
483 greg 2.1 rayorigin(&rr, PRIMARY, NULL, NULL);
484     if (ray_pqueue(&rr) == 1)
485     addcolor(vsum, rr.rcol);
486     }
487     /* finish MC calculation */
488     while (ray_presult(&rr, 0) > 0)
489     addcolor(vsum, rr.rcol);
490     scalecolor(vsum, gscale/(nt*np));
491     /* compute direct component */
492     for (i = ndirs; i-- > 0; ) {
493     SRCINDEX si;
494     initsrcindex(&si);
495     while (srcray(&rr, NULL, &si)) {
496     double d = sens_val(rr.rdir);
497     if (d <= FTINY)
498     continue;
499     d *= si.dom/ndirs;
500     scalecolor(rr.rcoef, d);
501     if (ray_pqueue(&rr) == 1) {
502     multcolor(rr.rcol, rr.rcoef);
503     addcolor(vsum, rr.rcol);
504     }
505     }
506     }
507     /* finish direct calculation */
508     while (ray_presult(&rr, 0) > 0) {
509     multcolor(rr.rcol, rr.rcoef);
510     addcolor(vsum, rr.rcol);
511     }
512     /* print our result */
513     printf("%.4e %.4e %.4e\n", colval(vsum,RED),
514     colval(vsum,GRN), colval(vsum,BLU));
515     }