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 (16 years, 7 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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rsensor.c,v 2.4 2008/02/22 21:52:10 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 int
75 main(
76 int argc,
77 char *argv[]
78 )
79 {
80 int doheader = 1;
81 int optwarn = 0;
82 int i, rval;
83
84 progname = argv[0];
85 /* set up rendering defaults */
86 rand_samp = 1;
87 dstrsrc = 0.5;
88 srcsizerat = 0.1;
89 directrelay = 3;
90 ambounce = 1;
91 maxdepth = -10;
92 /* get options from command line */
93 for (i = 1; i < argc; i++) {
94 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 if (argv[i][0] != '-') {
101 if (i >= argc-1)
102 break; /* final octree argument */
103 if (!ray_pnprocs) {
104 over_options();
105 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 comp_sensor(argv[i]); /* process a sensor file */
114 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 if (!strcmp(argv[i], "-w")) { /* toggle warnings */
157 nowarn = !nowarn;
158 continue;
159 }
160 if (ray_pnprocs) {
161 if (!optwarn++)
162 error(WARNING,
163 "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 if (!ray_pnprocs)
186 error(USER, i<argc ? "missing sensor file" : "missing octree");
187 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 sensor = NULL;
299 sntp[0] = sntp[1] = 0;
300 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 maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
318 thdiv[0] = .0;
319 for (t = 1; t < sntp[0]; t++)
320 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
321 thdiv[sntp[0]] = maxtheta*DEGREE;
322 phdiv[0] = .0;
323 for (p = 1; p < sntp[1]; p++)
324 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
325 phdiv[sntp[1]] = 2.*PI;
326 /* size our table */
327 tsize = 1. - cos(maxtheta*DEGREE);
328 psize = PI*tsize/(maxtheta*DEGREE);
329 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
330 samptot = sntp[0]*sntp[1];
331 ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
332 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 rowp = &s_val(i,0);
341 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 for (i = 0; i < ntheta; i++) {
349 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 /* 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 pvals[i*(nphi+1)] = .0f;
366 for (j = 1; j < nphi; j++) {
367 prob = (double)j / (double)nphi;
368 rowp = &s_val(t,0);
369 rowp1 = &s_val(t+1,0);
370 for (p = 0; p < sntp[1]; p++) {
371 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
372 frac*rowp1[p]/rowsum[t+1]) <= .0)
373 break;
374 if (p >= sntp[1])
375 error(INTERNAL,
376 "code error 2 in init_ptable()");
377 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
378 + frac*rowp1[p]/rowsum[t+1]);
379 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 tvals[0] = .0f;
386 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 while (t+1 < sntp[0] && theta >= s_theta(t+1))
442 ++t;
443 while (t-1 >= 0 && theta <= s_theta(t-1))
444 --t;
445 while (p+1 < sntp[1] && phi >= s_phi(p+1))
446 ++p;
447 while (p-1 >= 0 && phi <= s_phi(p-1))
448 --p;
449 return(s_val(t,p));
450 }
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 for (j = 0; j < np; j++) {
480 VCOPY(rr.rorg, ourview.vp);
481 get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
482 if (ourview.vfore > FTINY)
483 VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
484 rr.rmax = .0;
485 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 }