ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.18
Committed: Tue Jun 9 21:34:15 2015 UTC (8 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2, rad5R0, rad5R1
Changes since 2.17: +8 -4 lines
Log Message:
Added checks for extra and missing angles and documented behavior

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rsensor.c,v 2.17 2015/04/24 19:17:12 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 181 /* 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 prob. values (1-D table of 1-cos(t)) */
38 float *pvals = NULL; /* phi prob. 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 if (strcmp(argv[argc-1], "."))
124 ray_pinit(argv[argc-1], nprocs);
125 }
126 comp_sensor(argv[i]); /* process a sensor file */
127 continue;
128 }
129 if (argv[i][1] == 'r') { /* sampling options */
130 if (argv[i][2] == 'd')
131 nsamps = atol(argv[++i]);
132 else if (argv[i][2] == 's')
133 nssamps = atol(argv[++i]);
134 else {
135 sprintf(errmsg, "bad option at '%s'", argv[i]);
136 error(USER, errmsg);
137 }
138 continue;
139 }
140 /* direct component samples */
141 if (argv[i][1] == 'd' && argv[i][2] == 'n') {
142 ndsamps = atoi(argv[++i]);
143 continue;
144 }
145 if (argv[i][1] == 'v') { /* next sensor view */
146 if (argv[i][2] == 'f') {
147 rval = viewfile(argv[++i], &ourview, NULL);
148 if (rval < 0) {
149 sprintf(errmsg,
150 "cannot open view file \"%s\"",
151 argv[i]);
152 error(SYSTEM, errmsg);
153 } else if (rval == 0) {
154 sprintf(errmsg,
155 "bad view file \"%s\"",
156 argv[i]);
157 error(USER, errmsg);
158 }
159 continue;
160 }
161 rval = getviewopt(&ourview, argc-i, argv+i);
162 if (rval >= 0) {
163 i += rval;
164 continue;
165 }
166 sprintf(errmsg, "bad view option at '%s'", argv[i]);
167 error(USER, errmsg);
168 }
169 if (!strcmp(argv[i], "-w")) { /* toggle warnings */
170 nowarn = !nowarn;
171 continue;
172 }
173 if (ray_pnprocs) {
174 if (!optwarn++)
175 error(WARNING,
176 "rendering options should appear before first sensor");
177 } else if (!strcmp(argv[i], "-defaults")) {
178 print_defaults();
179 return(0);
180 }
181 if (argv[i][1] == 'h') { /* header toggle */
182 doheader = !doheader;
183 continue;
184 }
185 if (!strcmp(argv[i], "-n")) { /* number of processes */
186 nprocs = atoi(argv[++i]);
187 if (nprocs <= 0)
188 error(USER, "illegal number of processes");
189 continue;
190 }
191 rval = getrenderopt(argc-i, argv+i);
192 if (rval < 0) {
193 sprintf(errmsg, "bad render option at '%s'", argv[i]);
194 error(USER, errmsg);
195 }
196 i += rval;
197 }
198 if (sensor == NULL)
199 error(USER, i<argc ? "missing sensor file" : "missing octree");
200 quit(0);
201 }
202
203 /* Load sensor sensitivities (first row and column are angles) */
204 static float *
205 load_sensor(
206 int ntp[2],
207 char *sfile
208 )
209 {
210 int warnedneg;
211 char linebuf[8192];
212 int last_pos_val = 0;
213 int nelem = 1000;
214 float *sarr = (float *)malloc(sizeof(float)*nelem);
215 FILE *fp;
216 char *cp;
217 int i;
218
219 fp = frlibopen(sfile);
220 if (fp == NULL) {
221 sprintf(errmsg, "cannot open sensor file '%s'", sfile);
222 error(SYSTEM, errmsg);
223 }
224 fgets(linebuf, sizeof(linebuf), fp);
225 if (!strncmp(linebuf, "Elevation ", 10))
226 fgets(linebuf, sizeof(linebuf), fp);
227 /* get phi values */
228 sarr[0] = .0f;
229 if (strncmp(linebuf, "degrees", 7)) {
230 sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
231 error(USER, errmsg);
232 }
233 cp = sskip(linebuf);
234 ntp[1] = 0;
235 for ( ; ; ) {
236 sarr[ntp[1]+1] = atof(cp);
237 cp = fskip(cp);
238 if (cp == NULL)
239 break;
240 if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]+FTINY) {
241 sprintf(errmsg,
242 "Phi values not monotinically increasing in sensor file '%s'",
243 sfile);
244 error(USER, errmsg);
245 }
246 ++ntp[1];
247 }
248 warnedneg = 0;
249 ntp[0] = 0; /* get thetas + data */
250 while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
251 ++ntp[0];
252 if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
253 nelem += (nelem>>2) + ntp[1];
254 sarr = (float *)realloc((void *)sarr,
255 sizeof(float)*nelem);
256 if (sarr == NULL)
257 error(SYSTEM, "out of memory in load_sensor()");
258 }
259 cp = linebuf;
260 i = ntp[0]*(ntp[1]+1);
261 for ( ; ; ) {
262 sarr[i] = atof(cp);
263 cp = fskip(cp);
264 if (cp == NULL)
265 break;
266 if (sarr[i] < .0) {
267 if (!warnedneg++) {
268 sprintf(errmsg,
269 "Negative value(s) in sensor file '%s' (ignored)\n", sfile);
270 error(WARNING, errmsg);
271 }
272 sarr[i] = .0;
273 } else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1))
274 last_pos_val = i;
275 ++i;
276 }
277 if (i == ntp[0]*(ntp[1]+1)) /* empty line? */
278 break;
279 if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
280 sarr[(ntp[0]-1)*(ntp[1]+1)]) {
281 sprintf(errmsg,
282 "Theta values not monotinically increasing in sensor file '%s'",
283 sfile);
284 error(USER, errmsg);
285 }
286 if (i != (ntp[0]+1)*(ntp[1]+1)) {
287 sprintf(errmsg,
288 "bad column count near line %d in sensor file '%s'",
289 ntp[0]+1, sfile);
290 error(USER, errmsg);
291 }
292 }
293 /* truncate zero region */
294 ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1;
295 nelem = (ntp[0]+1)*(ntp[1]+1);
296 fclose(fp);
297 errmsg[0] = '\0'; /* sanity checks */
298 if (!last_pos_val)
299 sprintf(errmsg, "no positive sensor values in file '%s'", sfile);
300 else if (fabs(sarr[ntp[1]+1]) > FTINY)
301 sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
302 sfile);
303 else if (fabs(sarr[1]) > FTINY)
304 sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
305 sfile);
306 else if (sarr[ntp[1]] < 270.-FTINY)
307 sprintf(errmsg,
308 "maximum phi must be 270 or greater in sensor file '%s'",
309 sfile);
310 else if (sarr[ntp[1]] >= 360.-FTINY)
311 sprintf(errmsg,
312 "maximum phi must be less than 360 in sensor file '%s'",
313 sfile);
314 if (errmsg[0])
315 error(USER, errmsg);
316 return((float *)realloc((void *)sarr, sizeof(float)*nelem));
317 }
318
319 /* Initialize probability table */
320 static void
321 init_ptable(
322 char *sfile
323 )
324 {
325 int samptot = nsamps;
326 float *rowp, *rowp1;
327 double rowsum[MAXNT], rowomega[MAXNT];
328 double thdiv[MAXNT+1], phdiv[MAXNP+1];
329 double tsize, psize;
330 double prob, frac, frac1;
331 int i, j, t, p;
332 /* free old table */
333 if (sensor != NULL)
334 free((void *)sensor);
335 if (pvals != NULL)
336 free((void *)pvals);
337 if (sfile == NULL || !*sfile) {
338 sensor = NULL;
339 sntp[0] = sntp[1] = 0;
340 pvals = NULL;
341 ntheta = nphi = 0;
342 return;
343 }
344 /* load sensor table */
345 sensor = load_sensor(sntp, sfile);
346 if (sntp[0] > MAXNT) {
347 sprintf(errmsg, "Too many theta rows in sensor file '%s'",
348 sfile);
349 error(INTERNAL, errmsg);
350 }
351 if (sntp[1] > MAXNP) {
352 sprintf(errmsg, "Too many phi columns in sensor file '%s'",
353 sfile);
354 error(INTERNAL, errmsg);
355 }
356 /* compute boundary angles */
357 maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
358 if (maxtheta > PI)
359 maxtheta = PI;
360 thdiv[0] = .0;
361 for (t = 1; t < sntp[0]; t++)
362 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
363 thdiv[sntp[0]] = maxtheta;
364 phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
365 for (p = 1; p < sntp[1]; p++)
366 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
367 phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
368 /* size our table */
369 tsize = 1. - cos(maxtheta);
370 psize = PI*tsize/maxtheta;
371 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
372 samptot = sntp[0]*sntp[1];
373 ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1;
374 if (ntheta > MAXNT)
375 ntheta = MAXNT;
376 nphi = samptot/ntheta;
377 pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
378 if (pvals == NULL)
379 error(SYSTEM, "out of memory in init_ptable()");
380 gscale = .0; /* compute our inverse table */
381 for (i = 0; i < sntp[0]; i++) {
382 rowp = &s_val(i,0);
383 rowsum[i] = 1e-20;
384 for (j = 0; j < sntp[1]; j++)
385 rowsum[i] += *rowp++;
386 rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
387 rowomega[i] *= 2.*PI / (double)sntp[1];
388 gscale += rowsum[i] * rowomega[i];
389 }
390 if (gscale <= FTINY) {
391 sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
392 error(USER, errmsg);
393 }
394 for (i = 0; i < ntheta; i++) {
395 prob = (double)i / (double)ntheta;
396 for (t = 0; t < sntp[0]; t++)
397 if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
398 break;
399 if (t >= sntp[0])
400 error(INTERNAL, "code error 1 in init_ptable()");
401 frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
402 tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
403 frac*cos(thdiv[t+1]) );
404 /* offset b/c sensor values are centered */
405 if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) {
406 frac -= 0.5;
407 } else {
408 frac += 0.5;
409 --t;
410 }
411 pvals[i*(nphi+1)] = phdiv[0];
412 for (j = 1; j < nphi; j++) {
413 prob = (double)j / (double)nphi;
414 rowp = &s_val(t,0);
415 rowp1 = &s_val(t+1,0);
416 for (p = 0; p < sntp[1]; p++)
417 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
418 frac*rowp1[p]/rowsum[t+1]) <= .0)
419 break;
420 if (p >= sntp[1]) { /* should never happen? */
421 p = sntp[1] - 1;
422 prob = .5;
423 }
424 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
425 + frac*rowp1[p]/rowsum[t+1]);
426 pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
427 frac1*phdiv[p+1];
428 }
429 pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
430 }
431 /* duplicate final row */
432 memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1),
433 sizeof(*pvals)*(nphi+1));
434 tvals[0] = .0f;
435 tvals[ntheta] = (float)tsize;
436 }
437
438 /* Get normalized direction from random variables in [0,1) range */
439 static void
440 get_direc(
441 FVECT dvec,
442 double x,
443 double y
444 )
445 {
446 double xfrac = x*ntheta;
447 int tndx = (int)xfrac;
448 double yfrac = y*nphi;
449 int pndx = (int)yfrac;
450 double rad, phi;
451 FVECT dv;
452 int i;
453
454 xfrac -= (double)tndx;
455 yfrac -= (double)pndx;
456 pndx += tndx*(nphi+1);
457
458 dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]);
459 rad = sqrt(1. - dv[2]*dv[2]);
460 phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1];
461 dv[0] = -rad*sin(phi);
462 dv[1] = rad*cos(phi);
463 for (i = 3; i--; )
464 dvec[i] = dv[0]*ourview.hvec[i] +
465 dv[1]*ourview.vvec[i] +
466 dv[2]*ourview.vdir[i] ;
467 }
468
469 /* Get sensor value in the specified direction (normalized) */
470 static float
471 sens_val(
472 FVECT dvec
473 )
474 {
475 FVECT dv;
476 float theta, phi;
477 int t, p;
478
479 dv[2] = DOT(dvec, ourview.vdir);
480 theta = acos(dv[2]);
481 if (theta >= maxtheta)
482 return(.0f);
483 dv[0] = DOT(dvec, ourview.hvec);
484 dv[1] = DOT(dvec, ourview.vvec);
485 phi = atan2(-dv[0], dv[1]);
486 while (phi < .0f) phi += (float)(2.*PI);
487 t = (int)(theta/maxtheta * sntp[0]);
488 p = (int)(phi*(1./(2.*PI)) * sntp[1]);
489 /* hack for non-uniform sensor grid */
490 theta *= (float)(1./DEGREE);
491 phi *= (float)(1./DEGREE);
492 while (t+1 < sntp[0] && theta >= s_theta(t+1))
493 ++t;
494 while (t-1 >= 0 && theta <= s_theta(t-1))
495 --t;
496 while (p+1 < sntp[1] && phi >= s_phi(p+1))
497 ++p;
498 while (p-1 >= 0 && phi <= s_phi(p-1))
499 --p;
500 return(s_val(t,p));
501 }
502
503 /* Print origin and direction */
504 static void
505 print_ray(
506 FVECT rorg,
507 FVECT rdir
508 )
509 {
510 printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
511 rorg[0], rorg[1], rorg[2],
512 rdir[0], rdir[1], rdir[2]);
513 }
514
515 /* Compute sensor output */
516 static void
517 comp_sensor(
518 char *sfile
519 )
520 {
521 int ndirs = dstrsrc > FTINY ? ndsamps :
522 ndsamps > 0 ? 1 : 0;
523 char *err;
524 int nt, np;
525 COLOR vsum;
526 RAY rr;
527 double sf;
528 int i, j;
529 /* set view */
530 ourview.type = VT_ANG;
531 ourview.horiz = ourview.vert = 180.;
532 ourview.hoff = ourview.voff = .0;
533 err = setview(&ourview);
534 if (err != NULL)
535 error(USER, err);
536 /* assign probability table */
537 init_ptable(sfile);
538 /* stratified MC sampling */
539 setcolor(vsum, .0f, .0f, .0f);
540 nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
541 np = nsamps/nt;
542 sf = gscale/nsamps;
543 for (i = 0; i < nt; i++)
544 for (j = 0; j < np; j++) {
545 VCOPY(rr.rorg, ourview.vp);
546 get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
547 if (ourview.vfore > FTINY)
548 VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
549 if (!ray_pnprocs) {
550 print_ray(rr.rorg, rr.rdir);
551 continue;
552 }
553 rr.rmax = .0;
554 rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
555 scalecolor(rr.rcoef, sf);
556 if (ray_pqueue(&rr) == 1)
557 addcolor(vsum, rr.rcol);
558 }
559 /* remaining rays pure MC */
560 for (i = nsamps - nt*np; i-- > 0; ) {
561 VCOPY(rr.rorg, ourview.vp);
562 get_direc(rr.rdir, frandom(), frandom());
563 if (ourview.vfore > FTINY)
564 VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
565 if (!ray_pnprocs) {
566 print_ray(rr.rorg, rr.rdir);
567 continue;
568 }
569 rr.rmax = .0;
570 rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
571 scalecolor(rr.rcoef, sf);
572 if (ray_pqueue(&rr) == 1)
573 addcolor(vsum, rr.rcol);
574 }
575 if (!ray_pnprocs) /* just printing rays */
576 return;
577 /* scale partial result */
578 scalecolor(vsum, sf);
579 /* add direct component */
580 for (i = ndirs; i-- > 0; ) {
581 SRCINDEX si;
582 initsrcindex(&si);
583 while (srcray(&rr, NULL, &si)) {
584 sf = sens_val(rr.rdir);
585 if (sf <= FTINY)
586 continue;
587 sf *= si.dom/ndirs;
588 scalecolor(rr.rcoef, sf);
589 if (ray_pqueue(&rr) == 1) {
590 multcolor(rr.rcol, rr.rcoef);
591 addcolor(vsum, rr.rcol);
592 }
593 }
594 }
595 /* finish our calculation */
596 while (ray_presult(&rr, 0) > 0) {
597 multcolor(rr.rcol, rr.rcoef);
598 addcolor(vsum, rr.rcol);
599 }
600 /* print our result */
601 printf("%.4e %.4e %.4e\n", colval(vsum,RED),
602 colval(vsum,GRN), colval(vsum,BLU));
603 }