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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rsensor.c,v 2.6 2008/12/13 00:44:05 greg Exp $";
3 #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 = {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
29 unsigned long nsamps = 10000; /* desired number of initial samples */
30 unsigned long nssamps = 9000; /* number of super-samples */
31 int ndsamps = 32; /* number of direct samples */
32 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 #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 static void comp_sensor(char *sfile);
48
49 static void
50 over_options() /* overriding options */
51 {
52 directvis = (ndsamps <= 0);
53 do_irrad = 0;
54 }
55
56 static void
57 print_defaults() /* print out default parameters */
58 {
59 over_options();
60 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
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 int
86 main(
87 int argc,
88 char *argv[]
89 )
90 {
91 int doheader = 1;
92 int optwarn = 0;
93 int i, rval;
94
95 progname = argv[0];
96 /* set up rendering defaults */
97 rand_samp = 1;
98 dstrsrc = 0.65;
99 srcsizerat = 0.1;
100 directrelay = 3;
101 ambounce = 1;
102 maxdepth = -10;
103 /* get options from command line */
104 for (i = 1; i < argc; i++) {
105 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 if (argv[i][0] != '-') {
112 if (i >= argc-1)
113 break; /* final octree argument */
114 if (!ray_pnprocs) {
115 over_options();
116 if (doheader) { /* print header */
117 newheader("RADIANCE", stdout);
118 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 comp_sensor(argv[i]); /* process a sensor file */
126 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 if (!strcmp(argv[i], "-w")) { /* toggle warnings */
169 nowarn = !nowarn;
170 continue;
171 }
172 if (ray_pnprocs) {
173 if (!optwarn++)
174 error(WARNING,
175 "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 if (!ray_pnprocs)
198 error(USER, i<argc ? "missing sensor file" : "missing octree");
199 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 sensor = NULL;
311 sntp[0] = sntp[1] = 0;
312 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 maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
330 thdiv[0] = .0;
331 for (t = 1; t < sntp[0]; t++)
332 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
333 thdiv[sntp[0]] = maxtheta*DEGREE;
334 phdiv[0] = .0;
335 for (p = 1; p < sntp[1]; p++)
336 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
337 phdiv[sntp[1]] = 2.*PI;
338 /* size our table */
339 tsize = 1. - cos(maxtheta*DEGREE);
340 psize = PI*tsize/(maxtheta*DEGREE);
341 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
342 samptot = sntp[0]*sntp[1];
343 ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
344 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 rowp = &s_val(i,0);
353 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 for (i = 0; i < ntheta; i++) {
361 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 /* 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 pvals[i*(nphi+1)] = .0f;
378 for (j = 1; j < nphi; j++) {
379 prob = (double)j / (double)nphi;
380 rowp = &s_val(t,0);
381 rowp1 = &s_val(t+1,0);
382 for (p = 0; p < sntp[1]; p++) {
383 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
384 frac*rowp1[p]/rowsum[t+1]) <= .0)
385 break;
386 if (p >= sntp[1])
387 error(INTERNAL,
388 "code error 2 in init_ptable()");
389 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
390 + frac*rowp1[p]/rowsum[t+1]);
391 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 tvals[0] = .0f;
398 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 while (t+1 < sntp[0] && theta >= s_theta(t+1))
454 ++t;
455 while (t-1 >= 0 && theta <= s_theta(t-1))
456 --t;
457 while (p+1 < sntp[1] && phi >= s_phi(p+1))
458 ++p;
459 while (p-1 >= 0 && phi <= s_phi(p-1))
460 --p;
461 return(s_val(t,p));
462 }
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 double sf;
477 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 /* stratified MC sampling */
488 setcolor(vsum, .0f, .0f, .0f);
489 nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
490 np = nsamps/nt;
491 sf = gscale/nsamps;
492 for (i = 0; i < nt; i++)
493 for (j = 0; j < np; j++) {
494 VCOPY(rr.rorg, ourview.vp);
495 get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
496 if (ourview.vfore > FTINY)
497 VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
498 rr.rmax = .0;
499 rayorigin(&rr, PRIMARY, NULL, NULL);
500 scalecolor(rr.rcoef, sf);
501 if (ray_pqueue(&rr) == 1)
502 addcolor(vsum, rr.rcol);
503 }
504 /* 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 for (i = ndirs; i-- > 0; ) {
520 SRCINDEX si;
521 initsrcindex(&si);
522 while (srcray(&rr, NULL, &si)) {
523 sf = sens_val(rr.rdir);
524 if (sf <= FTINY)
525 continue;
526 sf *= si.dom/ndirs;
527 scalecolor(rr.rcoef, sf);
528 if (ray_pqueue(&rr) == 1) {
529 multcolor(rr.rcol, rr.rcoef);
530 addcolor(vsum, rr.rcol);
531 }
532 }
533 }
534 /* finish our calculation */
535 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 }