ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.21
Committed: Tue Mar 15 18:05:03 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4
Changes since 2.20: +3 -7 lines
Log Message:
fix(rsensor): Minor corrections wrt. integer types

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rsensor.c,v 2.20 2020/07/24 17:09:33 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 "platform.h"
13 #include "source.h"
14 #include "view.h"
15 #include "random.h"
16
17 #define DEGREE (PI/180.)
18
19 #define MAXNT 181 /* maximum number of theta divisions */
20 #define MAXNP 360 /* maximum number of phi divisions */
21
22 extern char *progname; /* global argv[0] */
23 extern int nowarn; /* don't report warnings? */
24
25 /* current sensor's perspective */
26 VIEW ourview = {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,0.,0.},
27 1.,180.,180.,0.,0.,0.,0.,
28 {0.,0.,0.},{0.,0.,0.},0.,0.};
29
30 long nsamps = 10000; /* desired number of initial 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("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
63 printf("-vp %f %f %f\t# view point\n",
64 ourview.vp[0], ourview.vp[1], ourview.vp[2]);
65 printf("-vd %f %f %f\t# view direction\n",
66 ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]);
67 printf("-vu %f %f %f\t# view up\n",
68 ourview.vup[0], ourview.vup[1], ourview.vup[2]);
69 printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore);
70 print_rdefaults();
71 }
72
73
74 void
75 quit(ec) /* make sure exit is called */
76 int ec;
77 {
78 if (ray_pnprocs > 0) /* close children if any */
79 ray_pclose(0);
80 else if (ray_pnprocs < 0)
81 _exit(ec); /* avoid flush in child */
82 exit(ec);
83 }
84
85
86 int
87 main(
88 int argc,
89 char *argv[]
90 )
91 {
92 int doheader = 1;
93 int optwarn = 0;
94 int i, rval;
95
96 progname = argv[0];
97 /* set up rendering defaults */
98 rand_samp = 1;
99 dstrsrc = 0.65;
100 srcsizerat = 0.1;
101 directrelay = 3;
102 ambounce = 1;
103 maxdepth = -10;
104 /* get options from command line */
105 for (i = 1; i < argc; i++) {
106 while ((rval = expandarg(&argc, &argv, i)) > 0)
107 ;
108 if (rval < 0) {
109 sprintf(errmsg, "cannot expand '%s'", argv[i]);
110 error(SYSTEM, errmsg);
111 }
112 if (argv[i][0] != '-') {
113 if (i >= argc-1)
114 break; /* final octree argument */
115 if (!ray_pnprocs) {
116 over_options();
117 if (doheader) { /* print header */
118 newheader("RADIANCE", stdout);
119 printargs(argc, argv, stdout);
120 fputformat("ascii", stdout);
121 putchar('\n');
122 }
123 /* start process(es) */
124 if (strcmp(argv[argc-1], "."))
125 ray_pinit(argv[argc-1], nprocs);
126 }
127 comp_sensor(argv[i]); /* process a sensor file */
128 continue;
129 }
130 if (argv[i][1] == 'r') { /* sampling options */
131 if (argv[i][2] == 'd')
132 nsamps = 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 (sensor == NULL)
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 int warnedneg;
210 char linebuf[8192];
211 int last_pos_val = 0;
212 int nelem = 1000;
213 float *sarr = (float *)malloc(sizeof(float)*nelem);
214 FILE *fp;
215 char *cp;
216 int i;
217
218 fp = frlibopen(sfile);
219 if (fp == NULL) {
220 sprintf(errmsg, "cannot open sensor file '%s'", sfile);
221 error(SYSTEM, errmsg);
222 }
223 fgets(linebuf, sizeof(linebuf), fp);
224 if (!strncmp(linebuf, "Elevation ", 10))
225 fgets(linebuf, sizeof(linebuf), fp);
226 /* get phi values */
227 sarr[0] = .0f;
228 if (strncmp(linebuf, "degrees", 7)) {
229 sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
230 error(USER, errmsg);
231 }
232 cp = sskip(linebuf);
233 ntp[1] = 0;
234 for ( ; ; ) {
235 sarr[ntp[1]+1] = atof(cp);
236 cp = fskip(cp);
237 if (cp == NULL)
238 break;
239 if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]+FTINY) {
240 sprintf(errmsg,
241 "Phi values not monotinically increasing in sensor file '%s'",
242 sfile);
243 error(USER, errmsg);
244 }
245 ++ntp[1];
246 }
247 warnedneg = 0;
248 ntp[0] = 0; /* get thetas + data */
249 while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
250 ++ntp[0];
251 if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
252 nelem += (nelem>>2) + ntp[1];
253 sarr = (float *)realloc((void *)sarr,
254 sizeof(float)*nelem);
255 if (sarr == NULL)
256 error(SYSTEM, "out of memory in load_sensor()");
257 }
258 cp = linebuf;
259 i = ntp[0]*(ntp[1]+1);
260 for ( ; ; ) {
261 sarr[i] = atof(cp);
262 cp = fskip(cp);
263 if (cp == NULL)
264 break;
265 if (sarr[i] < .0) {
266 if (!warnedneg++) {
267 sprintf(errmsg,
268 "Negative value(s) in sensor file '%s' (ignored)\n", sfile);
269 error(WARNING, errmsg);
270 }
271 sarr[i] = .0;
272 } else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1))
273 last_pos_val = i;
274 ++i;
275 }
276 if (i == ntp[0]*(ntp[1]+1)) /* empty line? */
277 break;
278 if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
279 sarr[(ntp[0]-1)*(ntp[1]+1)]) {
280 sprintf(errmsg,
281 "Theta values not monotinically increasing in sensor file '%s'",
282 sfile);
283 error(USER, errmsg);
284 }
285 if (i != (ntp[0]+1)*(ntp[1]+1)) {
286 sprintf(errmsg,
287 "bad column count near line %d in sensor file '%s'",
288 ntp[0]+1, sfile);
289 error(USER, errmsg);
290 }
291 }
292 /* truncate zero region */
293 ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1;
294 nelem = (ntp[0]+1)*(ntp[1]+1);
295 fclose(fp);
296 errmsg[0] = '\0'; /* sanity checks */
297 if (!last_pos_val)
298 sprintf(errmsg, "no positive sensor values in file '%s'", sfile);
299 else if (fabs(sarr[ntp[1]+1]) > FTINY)
300 sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
301 sfile);
302 else if (fabs(sarr[1]) > FTINY)
303 sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
304 sfile);
305 else if (sarr[ntp[1]] < 270.-FTINY)
306 sprintf(errmsg,
307 "maximum phi must be 270 or greater in sensor file '%s'",
308 sfile);
309 else if (sarr[ntp[1]] >= 360.-FTINY)
310 sprintf(errmsg,
311 "maximum phi must be less than 360 in sensor file '%s'",
312 sfile);
313 if (errmsg[0])
314 error(USER, errmsg);
315 return((float *)realloc((void *)sarr, sizeof(float)*nelem));
316 }
317
318 /* Initialize probability table */
319 static void
320 init_ptable(
321 char *sfile
322 )
323 {
324 long samptot = nsamps;
325 float *rowp, *rowp1;
326 double rowsum[MAXNT], rowomega[MAXNT];
327 double thdiv[MAXNT+1], phdiv[MAXNP+1];
328 double tsize, psize;
329 double prob, frac, frac1;
330 int i, j, t, p;
331 /* free old table */
332 if (sensor != NULL)
333 free((void *)sensor);
334 if (pvals != NULL)
335 free((void *)pvals);
336 if (sfile == NULL || !*sfile) {
337 sensor = NULL;
338 sntp[0] = sntp[1] = 0;
339 pvals = NULL;
340 ntheta = nphi = 0;
341 return;
342 }
343 /* load sensor table */
344 sensor = load_sensor(sntp, sfile);
345 if (sntp[0] > MAXNT) {
346 sprintf(errmsg, "Too many theta rows in sensor file '%s'",
347 sfile);
348 error(INTERNAL, errmsg);
349 }
350 if (sntp[1] > MAXNP) {
351 sprintf(errmsg, "Too many phi columns in sensor file '%s'",
352 sfile);
353 error(INTERNAL, errmsg);
354 }
355 /* compute boundary angles */
356 maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
357 if (maxtheta > PI)
358 maxtheta = PI;
359 thdiv[0] = .0;
360 for (t = 1; t < sntp[0]; t++)
361 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
362 thdiv[sntp[0]] = maxtheta;
363 phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
364 for (p = 1; p < sntp[1]; p++)
365 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
366 phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
367 /* size our table */
368 tsize = 1. - cos(maxtheta);
369 psize = PI*tsize/maxtheta;
370 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
371 samptot = sntp[0]*sntp[1];
372 ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1;
373 if (ntheta > MAXNT)
374 ntheta = MAXNT;
375 nphi = samptot/ntheta;
376 pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
377 if (pvals == NULL)
378 error(SYSTEM, "out of memory in init_ptable()");
379 gscale = .0; /* compute our inverse table */
380 for (i = 0; i < sntp[0]; i++) {
381 rowp = &s_val(i,0);
382 rowsum[i] = 1e-20;
383 for (j = 0; j < sntp[1]; j++)
384 rowsum[i] += *rowp++;
385 rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
386 rowomega[i] *= 2.*PI / (double)sntp[1];
387 gscale += rowsum[i] * rowomega[i];
388 }
389 if (gscale <= FTINY) {
390 sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
391 error(USER, errmsg);
392 }
393 for (i = 0; i < ntheta; i++) {
394 prob = (double)i / (double)ntheta;
395 for (t = 0; t < sntp[0]; t++)
396 if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
397 break;
398 if (t >= sntp[0])
399 error(INTERNAL, "code error 1 in init_ptable()");
400 frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
401 tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
402 frac*cos(thdiv[t+1]) );
403 /* offset b/c sensor values are centered */
404 if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) {
405 frac -= 0.5;
406 } else {
407 frac += 0.5;
408 --t;
409 }
410 pvals[i*(nphi+1)] = phdiv[0];
411 for (j = 1; j < nphi; j++) {
412 prob = (double)j / (double)nphi;
413 rowp = &s_val(t,0);
414 rowp1 = &s_val(t+1,0);
415 for (p = 0; p < sntp[1]; p++)
416 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
417 frac*rowp1[p]/rowsum[t+1]) <= .0)
418 break;
419 if (p >= sntp[1]) { /* should never happen? */
420 p = sntp[1] - 1;
421 prob = .5;
422 }
423 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
424 + frac*rowp1[p]/rowsum[t+1]);
425 pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
426 frac1*phdiv[p+1];
427 }
428 pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
429 }
430 /* duplicate final row */
431 memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1),
432 sizeof(*pvals)*(nphi+1));
433 tvals[0] = .0f;
434 tvals[ntheta] = (float)tsize;
435 }
436
437 /* Get normalized direction from random variables in [0,1) range */
438 static void
439 get_direc(
440 FVECT dvec,
441 double x,
442 double y
443 )
444 {
445 double xfrac = x*ntheta;
446 int tndx = (int)xfrac;
447 double yfrac = y*nphi;
448 int pndx = (int)yfrac;
449 double rad, phi;
450 FVECT dv;
451 int i;
452
453 xfrac -= (double)tndx;
454 yfrac -= (double)pndx;
455 pndx += tndx*(nphi+1);
456
457 dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]);
458 rad = sqrt(1. - dv[2]*dv[2]);
459 phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1];
460 dv[0] = -rad*sin(phi);
461 dv[1] = rad*cos(phi);
462 for (i = 3; i--; )
463 dvec[i] = dv[0]*ourview.hvec[i] +
464 dv[1]*ourview.vvec[i] +
465 dv[2]*ourview.vdir[i] ;
466 }
467
468 /* Get sensor value in the specified direction (normalized) */
469 static float
470 sens_val(
471 FVECT dvec
472 )
473 {
474 FVECT dv;
475 float theta, phi;
476 int t, p;
477
478 dv[2] = DOT(dvec, ourview.vdir);
479 theta = acos(dv[2]);
480 if (theta >= maxtheta)
481 return(.0f);
482 dv[0] = DOT(dvec, ourview.hvec);
483 dv[1] = DOT(dvec, ourview.vvec);
484 phi = atan2(-dv[0], dv[1]);
485 while (phi < .0f) phi += (float)(2.*PI);
486 t = (int)(theta/maxtheta * sntp[0]);
487 p = (int)(phi*(1./(2.*PI)) * sntp[1]);
488 /* hack for non-uniform sensor grid */
489 theta *= (float)(1./DEGREE);
490 phi *= (float)(1./DEGREE);
491 while (t+1 < sntp[0] && theta >= s_theta(t+1))
492 ++t;
493 while (t-1 >= 0 && theta <= s_theta(t-1))
494 --t;
495 while (p+1 < sntp[1] && phi >= s_phi(p+1))
496 ++p;
497 while (p-1 >= 0 && phi <= s_phi(p-1))
498 --p;
499 return(s_val(t,p));
500 }
501
502 /* Print origin and direction */
503 static void
504 print_ray(
505 FVECT rorg,
506 FVECT rdir
507 )
508 {
509 printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
510 rorg[0], rorg[1], rorg[2],
511 rdir[0], rdir[1], rdir[2]);
512 }
513
514 /* Compute sensor output */
515 static void
516 comp_sensor(
517 char *sfile
518 )
519 {
520 int ndirs = dstrsrc > FTINY ? ndsamps :
521 ndsamps > 0 ? 1 : 0;
522 char *err;
523 int nt, np;
524 COLOR vsum;
525 RAY rr;
526 double sf;
527 int i, j;
528 /* set view */
529 ourview.type = VT_ANG;
530 ourview.horiz = ourview.vert = 180.;
531 ourview.hoff = ourview.voff = .0;
532 err = setview(&ourview);
533 if (err != NULL)
534 error(USER, err);
535 /* assign probability table */
536 init_ptable(sfile);
537 /* stratified MC sampling */
538 setcolor(vsum, .0f, .0f, .0f);
539 nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
540 np = nsamps/nt;
541 sf = gscale/nsamps;
542 for (i = 0; i < nt; i++)
543 for (j = 0; j < np; j++) {
544 VCOPY(rr.rorg, ourview.vp);
545 get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
546 if (ourview.vfore > FTINY)
547 VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
548 if (!ray_pnprocs) {
549 print_ray(rr.rorg, rr.rdir);
550 continue;
551 }
552 rr.rmax = .0;
553 rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
554 scalecolor(rr.rcoef, sf);
555 if (ray_pqueue(&rr) == 1)
556 addcolor(vsum, rr.rcol);
557 }
558 /* remaining rays pure MC */
559 for (i = nsamps - nt*np; i-- > 0; ) {
560 VCOPY(rr.rorg, ourview.vp);
561 get_direc(rr.rdir, frandom(), frandom());
562 if (ourview.vfore > FTINY)
563 VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
564 if (!ray_pnprocs) {
565 print_ray(rr.rorg, rr.rdir);
566 continue;
567 }
568 rr.rmax = .0;
569 rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
570 scalecolor(rr.rcoef, sf);
571 if (ray_pqueue(&rr) == 1)
572 addcolor(vsum, rr.rcol);
573 }
574 if (!ray_pnprocs) /* just printing rays */
575 return;
576 /* scale partial result */
577 scalecolor(vsum, sf);
578 /* add direct component */
579 for (i = ndirs; i-- > 0; ) {
580 SRCINDEX si;
581 initsrcindex(&si);
582 while (srcray(&rr, NULL, &si)) {
583 sf = sens_val(rr.rdir);
584 if (sf <= FTINY)
585 continue;
586 sf *= si.dom/ndirs;
587 scalecolor(rr.rcoef, sf);
588 if (ray_pqueue(&rr) == 1) {
589 multcolor(rr.rcol, rr.rcoef);
590 addcolor(vsum, rr.rcol);
591 }
592 }
593 }
594 /* finish our calculation */
595 while (ray_presult(&rr, 0) > 0) {
596 multcolor(rr.rcol, rr.rcoef);
597 addcolor(vsum, rr.rcol);
598 }
599 /* print our result */
600 printf("%.4e %.4e %.4e\n", colval(vsum,RED),
601 colval(vsum,GRN), colval(vsum,BLU));
602 }