ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.4
Committed: Fri Feb 22 21:52:10 2008 UTC (16 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.3: +1 -7 lines
Log Message:
Fixed minor inaccuracy in phi lookup

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rsensor.c,v 2.3 2008/02/22 18:04:02 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 = STDVIEW;
26
27 unsigned long nsamps = 10000; /* desired number of initial samples */
28 unsigned long nssamps = 9000; /* number of super-samples */
29 int ndsamps = 32; /* number of direct samples */
30 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 #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 static void comp_sensor(char *sfile);
46
47 static void
48 over_options() /* overriding options */
49 {
50 directvis = (ndsamps <= 0);
51 do_irrad = 0;
52 }
53
54 static void
55 print_defaults() /* print out default parameters */
56 {
57 over_options();
58 printf("-n %-9d\t\t\t# number of processes\n", nprocs);
59 printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
60 /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
61 printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
62 printf("-vp %f %f %f\t# view point\n",
63 ourview.vp[0], ourview.vp[1], ourview.vp[2]);
64 printf("-vd %f %f %f\t# view direction\n",
65 ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]);
66 printf("-vu %f %f %f\t# view up\n",
67 ourview.vup[0], ourview.vup[1], ourview.vup[2]);
68 printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore);
69 print_rdefaults();
70 }
71
72 int
73 main(
74 int argc,
75 char *argv[]
76 )
77 {
78 int doheader = 1;
79 int optwarn = 0;
80 int i, rval;
81
82 progname = argv[0];
83 /* set up rendering defaults */
84 rand_samp = 1;
85 dstrsrc = 0.5;
86 srcsizerat = 0.1;
87 directrelay = 3;
88 ambounce = 1;
89 maxdepth = -10;
90 /* get options from command line */
91 for (i = 1; i < argc; i++) {
92 while ((rval = expandarg(&argc, &argv, i)) > 0)
93 ;
94 if (rval < 0) {
95 sprintf(errmsg, "cannot expand '%s'", argv[i]);
96 error(SYSTEM, errmsg);
97 }
98 if (argv[i][0] != '-') {
99 if (i >= argc-1)
100 break; /* final octree argument */
101 if (!ray_pnprocs) {
102 over_options();
103 if (doheader) { /* print header */
104 printargs(argc, argv, stdout);
105 fputformat("ascii", stdout);
106 putchar('\n');
107 }
108 /* start process(es) */
109 ray_pinit(argv[argc-1], nprocs);
110 }
111 comp_sensor(argv[i]); /* process a sensor file */
112 continue;
113 }
114 if (argv[i][1] == 'r') { /* sampling options */
115 if (argv[i][2] == 'd')
116 nsamps = atol(argv[++i]);
117 else if (argv[i][2] == 's')
118 nssamps = atol(argv[++i]);
119 else {
120 sprintf(errmsg, "bad option at '%s'", argv[i]);
121 error(USER, errmsg);
122 }
123 continue;
124 }
125 /* direct component samples */
126 if (argv[i][1] == 'd' && argv[i][2] == 'n') {
127 ndsamps = atoi(argv[++i]);
128 continue;
129 }
130 if (argv[i][1] == 'v') { /* next sensor view */
131 if (argv[i][2] == 'f') {
132 rval = viewfile(argv[++i], &ourview, NULL);
133 if (rval < 0) {
134 sprintf(errmsg,
135 "cannot open view file \"%s\"",
136 argv[i]);
137 error(SYSTEM, errmsg);
138 } else if (rval == 0) {
139 sprintf(errmsg,
140 "bad view file \"%s\"",
141 argv[i]);
142 error(USER, errmsg);
143 }
144 continue;
145 }
146 rval = getviewopt(&ourview, argc-i, argv+i);
147 if (rval >= 0) {
148 i += rval;
149 continue;
150 }
151 sprintf(errmsg, "bad view option at '%s'", argv[i]);
152 error(USER, errmsg);
153 }
154 if (!strcmp(argv[i], "-w")) { /* toggle warnings */
155 nowarn = !nowarn;
156 continue;
157 }
158 if (ray_pnprocs) {
159 if (!optwarn++)
160 error(WARNING,
161 "rendering options should appear before first sensor");
162 } else if (!strcmp(argv[i], "-defaults")) {
163 print_defaults();
164 return(0);
165 }
166 if (argv[i][1] == 'h') { /* header toggle */
167 doheader = !doheader;
168 continue;
169 }
170 if (!strcmp(argv[i], "-n")) { /* number of processes */
171 nprocs = atoi(argv[++i]);
172 if (nprocs <= 0)
173 error(USER, "illegal number of processes");
174 continue;
175 }
176 rval = getrenderopt(argc-i, argv+i);
177 if (rval < 0) {
178 sprintf(errmsg, "bad render option at '%s'", argv[i]);
179 error(USER, errmsg);
180 }
181 i += rval;
182 }
183 if (!ray_pnprocs)
184 error(USER, i<argc ? "missing sensor file" : "missing octree");
185 quit(0);
186 }
187
188 /* Load sensor sensitivities (first row and column are angles) */
189 static float *
190 load_sensor(
191 int ntp[2],
192 char *sfile
193 )
194 {
195 char linebuf[8192];
196 int nelem = 1000;
197 float *sarr = (float *)malloc(sizeof(float)*nelem);
198 FILE *fp;
199 char *cp;
200 int i;
201
202 fp = frlibopen(sfile);
203 if (fp == NULL) {
204 sprintf(errmsg, "cannot open sensor file '%s'", sfile);
205 error(SYSTEM, errmsg);
206 }
207 fgets(linebuf, sizeof(linebuf), fp);
208 if (!strncmp(linebuf, "Elevation ", 10))
209 fgets(linebuf, sizeof(linebuf), fp);
210 /* get phi values */
211 sarr[0] = .0f;
212 if (strncmp(linebuf, "degrees", 7)) {
213 sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
214 error(USER, errmsg);
215 }
216 cp = sskip(linebuf);
217 ntp[1] = 0;
218 for ( ; ; ) {
219 sarr[ntp[1]+1] = atof(cp);
220 cp = fskip(cp);
221 if (cp == NULL)
222 break;
223 ++ntp[1];
224 }
225 ntp[0] = 0; /* get thetas + data */
226 while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
227 ++ntp[0];
228 if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
229 nelem += (nelem>>2) + ntp[1];
230 sarr = (float *)realloc((void *)sarr,
231 sizeof(float)*nelem);
232 if (sarr == NULL)
233 error(SYSTEM, "out of memory in load_sensor()");
234 }
235 cp = linebuf;
236 i = ntp[0]*(ntp[1]+1);
237 for ( ; ; ) {
238 sarr[i] = atof(cp);
239 cp = fskip(cp);
240 if (cp == NULL)
241 break;
242 ++i;
243 }
244 if (i == ntp[0]*(ntp[1]+1))
245 break;
246 if (i != (ntp[0]+1)*(ntp[1]+1)) {
247 sprintf(errmsg,
248 "bad column count near line %d in sensor file '%s'",
249 ntp[0]+1, sfile);
250 error(USER, errmsg);
251 }
252 }
253 nelem = i;
254 fclose(fp);
255 errmsg[0] = '\0'; /* sanity checks */
256 if (ntp[0] <= 0)
257 sprintf(errmsg, "no data in sensor file '%s'", sfile);
258 else if (fabs(sarr[ntp[1]+1]) > FTINY)
259 sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
260 sfile);
261 else if (fabs(sarr[1]) > FTINY)
262 sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
263 sfile);
264 else if (sarr[ntp[1]] <= FTINY)
265 sprintf(errmsg,
266 "maximum phi must be positive in sensor file '%s'",
267 sfile);
268 else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
269 sprintf(errmsg,
270 "maximum theta must be positive in sensor file '%s'",
271 sfile);
272 if (errmsg[0])
273 error(USER, errmsg);
274 return((float *)realloc((void *)sarr, sizeof(float)*nelem));
275 }
276
277 /* Initialize probability table */
278 static void
279 init_ptable(
280 char *sfile
281 )
282 {
283 int samptot = nsamps;
284 float *rowp, *rowp1;
285 double rowsum[MAXNT], rowomega[MAXNT];
286 double thdiv[MAXNT+1], phdiv[MAXNP+1];
287 double tsize, psize;
288 double prob, frac, frac1;
289 int i, j, t, p;
290 /* free old table */
291 if (sensor != NULL)
292 free((void *)sensor);
293 if (pvals != NULL)
294 free((void *)pvals);
295 if (sfile == NULL || !*sfile) {
296 sensor = NULL;
297 sntp[0] = sntp[1] = 0;
298 pvals = NULL;
299 ntheta = nphi = 0;
300 return;
301 }
302 /* load sensor table */
303 sensor = load_sensor(sntp, sfile);
304 if (sntp[0] > MAXNT) {
305 sprintf(errmsg, "Too many theta rows in sensor file '%s'",
306 sfile);
307 error(INTERNAL, errmsg);
308 }
309 if (sntp[1] > MAXNP) {
310 sprintf(errmsg, "Too many phi columns in sensor file '%s'",
311 sfile);
312 error(INTERNAL, errmsg);
313 }
314 /* compute boundary angles */
315 maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
316 thdiv[0] = .0;
317 for (t = 1; t < sntp[0]; t++)
318 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
319 thdiv[sntp[0]] = maxtheta*DEGREE;
320 phdiv[0] = .0;
321 for (p = 1; p < sntp[1]; p++)
322 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
323 phdiv[sntp[1]] = 2.*PI;
324 /* size our table */
325 tsize = 1. - cos(maxtheta*DEGREE);
326 psize = PI*tsize/(maxtheta*DEGREE);
327 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
328 samptot = sntp[0]*sntp[1];
329 ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
330 if (ntheta > MAXNT)
331 ntheta = MAXNT;
332 nphi = samptot/ntheta;
333 pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
334 if (pvals == NULL)
335 error(SYSTEM, "out of memory in init_ptable()");
336 gscale = .0; /* compute our inverse table */
337 for (i = 0; i < sntp[0]; i++) {
338 rowp = &s_val(i,0);
339 rowsum[i] = 0.;
340 for (j = 0; j < sntp[1]; j++)
341 rowsum[i] += *rowp++;
342 rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
343 rowomega[i] *= 2.*PI / (double)sntp[1];
344 gscale += rowsum[i] * rowomega[i];
345 }
346 for (i = 0; i < ntheta; i++) {
347 prob = (double)i / (double)ntheta;
348 for (t = 0; t < sntp[0]; t++)
349 if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
350 break;
351 if (t >= sntp[0])
352 error(INTERNAL, "code error 1 in init_ptable()");
353 frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
354 tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
355 frac*cos(thdiv[t+1]) );
356 /* offset b/c sensor values are centered */
357 if (t <= 0 || frac > 0.5)
358 frac -= 0.5;
359 else if (t >= sntp[0]-1 || frac < 0.5) {
360 frac += 0.5;
361 --t;
362 }
363 pvals[i*(nphi+1)] = .0f;
364 for (j = 1; j < nphi; j++) {
365 prob = (double)j / (double)nphi;
366 rowp = &s_val(t,0);
367 rowp1 = &s_val(t+1,0);
368 for (p = 0; p < sntp[1]; p++) {
369 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
370 frac*rowp1[p]/rowsum[t+1]) <= .0)
371 break;
372 if (p >= sntp[1])
373 error(INTERNAL,
374 "code error 2 in init_ptable()");
375 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
376 + frac*rowp1[p]/rowsum[t+1]);
377 pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
378 frac1*phdiv[p+1];
379 }
380 }
381 pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
382 }
383 tvals[0] = .0f;
384 tvals[ntheta] = (float)tsize;
385 }
386
387 /* Get normalized direction from random variables in [0,1) range */
388 static void
389 get_direc(
390 FVECT dvec,
391 double x,
392 double y
393 )
394 {
395 double xfrac = x*ntheta;
396 int tndx = (int)xfrac;
397 double yfrac = y*nphi;
398 int pndx = (int)yfrac;
399 double rad, phi;
400 FVECT dv;
401 int i;
402
403 xfrac -= (double)tndx;
404 yfrac -= (double)pndx;
405 pndx += tndx*(nphi+1);
406
407 dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]);
408 rad = sqrt(1. - dv[2]*dv[2]);
409 phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1];
410 dv[0] = -rad*sin(phi);
411 dv[1] = rad*cos(phi);
412 for (i = 3; i--; )
413 dvec[i] = dv[0]*ourview.hvec[i] +
414 dv[1]*ourview.vvec[i] +
415 dv[2]*ourview.vdir[i] ;
416 }
417
418 /* Get sensor value in the specified direction (normalized) */
419 static float
420 sens_val(
421 FVECT dvec
422 )
423 {
424 FVECT dv;
425 float theta, phi;
426 int t, p;
427
428 dv[2] = DOT(dvec, ourview.vdir);
429 theta = (float)((1./DEGREE) * acos(dv[2]));
430 if (theta >= maxtheta)
431 return(.0f);
432 dv[0] = DOT(dvec, ourview.hvec);
433 dv[1] = DOT(dvec, ourview.vvec);
434 phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
435 while (phi < .0f) phi += 360.f;
436 t = (int)(theta/maxtheta * sntp[0]);
437 p = (int)(phi*(1./360.) * sntp[1]);
438 /* hack for non-uniform sensor grid */
439 while (t+1 < sntp[0] && theta >= s_theta(t+1))
440 ++t;
441 while (t-1 >= 0 && theta <= s_theta(t-1))
442 --t;
443 while (p+1 < sntp[1] && phi >= s_phi(p+1))
444 ++p;
445 while (p-1 >= 0 && phi <= s_phi(p-1))
446 --p;
447 return(s_val(t,p));
448 }
449
450 /* Compute sensor output */
451 static void
452 comp_sensor(
453 char *sfile
454 )
455 {
456 int ndirs = dstrsrc > FTINY ? ndsamps :
457 ndsamps > 0 ? 1 : 0;
458 char *err;
459 int nt, np;
460 COLOR vsum;
461 RAY rr;
462 int i, j;
463 /* set view */
464 ourview.type = VT_ANG;
465 ourview.horiz = ourview.vert = 180.;
466 ourview.hoff = ourview.voff = .0;
467 err = setview(&ourview);
468 if (err != NULL)
469 error(USER, err);
470 /* assign probability table */
471 init_ptable(sfile);
472 /* do Monte Carlo sampling */
473 setcolor(vsum, .0f, .0f, .0f);
474 nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
475 np = nsamps/nt;
476 VCOPY(rr.rorg, ourview.vp);
477 rr.rmax = .0;
478 for (i = 0; i < nt; i++)
479 for (j = 0; j < np; j++) {
480 get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
481 rayorigin(&rr, PRIMARY, NULL, NULL);
482 if (ray_pqueue(&rr) == 1)
483 addcolor(vsum, rr.rcol);
484 }
485 /* finish MC calculation */
486 while (ray_presult(&rr, 0) > 0)
487 addcolor(vsum, rr.rcol);
488 scalecolor(vsum, gscale/(nt*np));
489 /* compute direct component */
490 for (i = ndirs; i-- > 0; ) {
491 SRCINDEX si;
492 initsrcindex(&si);
493 while (srcray(&rr, NULL, &si)) {
494 double d = sens_val(rr.rdir);
495 if (d <= FTINY)
496 continue;
497 d *= si.dom/ndirs;
498 scalecolor(rr.rcoef, d);
499 if (ray_pqueue(&rr) == 1) {
500 multcolor(rr.rcol, rr.rcoef);
501 addcolor(vsum, rr.rcol);
502 }
503 }
504 }
505 /* finish direct calculation */
506 while (ray_presult(&rr, 0) > 0) {
507 multcolor(rr.rcol, rr.rcoef);
508 addcolor(vsum, rr.rcol);
509 }
510 /* print our result */
511 printf("%.4e %.4e %.4e\n", colval(vsum,RED),
512 colval(vsum,GRN), colval(vsum,BLU));
513 }