ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.1
Committed: Thu Feb 21 01:22:06 2008 UTC (16 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial check-in of untested rsensor

File Contents

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