ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
(Generate patch)

Comparing ray/src/util/rsensor.c (file contents):
Revision 2.1 by greg, Thu Feb 21 01:22:06 2008 UTC vs.
Revision 2.9 by greg, Tue May 17 19:34:36 2011 UTC

# Line 22 | Line 22 | extern char    *progname;      /* global argv[0] */
22   extern int      nowarn;         /* don't report warnings? */
23  
24                                  /* current sensor's perspective */
25 < VIEW            ourview = STDVIEW;
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 = 16;   /* number of direct 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) */
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 < print_defaults()
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); */
# Line 57 | Line 71 | print_defaults()
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,
# Line 64 | Line 89 | main(
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 <        dstrsrc = 0.25;
97 >        rand_samp = 1;
98 >        dstrsrc = 0.65;
99 >        srcsizerat = 0.1;
100          directrelay = 3;
101          ambounce = 1;
102 <                                /* 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");
102 >        maxdepth = -10;
103                                  /* get options from command line */
104 <        for (i = 1; i < argc-1; i++) {
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] != '-') {        /* process a sensor file */
111 >                if (argv[i][0] != '-') {
112 >                        if (i >= argc-1)
113 >                                break;          /* final octree argument */
114                          if (!ray_pnprocs) {
115 <                                                /* overriding options */
93 <                                directvis = (ndsamps <= 0);
94 <                                do_irrad = 0;
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 <                                ray_pinit(argv[argc-1], nprocs);
123 >                                if (strcmp(argv[argc-1], "."))
124 >                                        ray_pinit(argv[argc-1], nprocs);
125                          }
126 <                        comp_sensor(argv[i]);
126 >                        comp_sensor(argv[i]);   /* process a sensor file */
127                          continue;
128                  }
129                  if (argv[i][1] == 'r') {        /* sampling options */
# Line 143 | Line 166 | main(
166                          sprintf(errmsg, "bad view option at '%s'", argv[i]);
167                          error(USER, errmsg);
168                  }
169 <                if (!strcmp(argv[i], "-w")) {   /* turn off warnings */
170 <                        nowarn = 1;
169 >                if (!strcmp(argv[i], "-w")) {   /* toggle warnings */
170 >                        nowarn = !nowarn;
171                          continue;
172                  }
173                  if (ray_pnprocs) {
174 <                        error(WARNING,
174 >                        if (!optwarn++)
175 >                                error(WARNING,
176                          "rendering options should appear before first sensor");
177                  } else if (!strcmp(argv[i], "-defaults")) {
178                          print_defaults();
# Line 171 | Line 195 | main(
195                  }
196                  i += rval;
197          }
198 +        if (sensor == NULL)
199 +                error(USER, i<argc ? "missing sensor file" : "missing octree");
200          quit(0);
201   }
202  
# Line 209 | Line 235 | load_sensor(
235                  cp = fskip(cp);
236                  if (cp == NULL)
237                          break;
238 +                if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]) {
239 +                        sprintf(errmsg,
240 +                "Phi values not monotinically increasing in sensor file '%s'",
241 +                                        sfile);
242 +                        error(USER, errmsg);
243 +                }
244                  ++ntp[1];
245          }
246          ntp[0] = 0;                             /* get thetas + data */
# Line 232 | Line 264 | load_sensor(
264                  }
265                  if (i == ntp[0]*(ntp[1]+1))
266                          break;
267 +                if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
268 +                                        sarr[(ntp[0]-1)*(ntp[1]+1)]) {
269 +                        sprintf(errmsg,
270 +                "Theta values not monotinically increasing in sensor file '%s'",
271 +                                        sfile);
272 +                        error(USER, errmsg);
273 +                }
274                  if (i != (ntp[0]+1)*(ntp[1]+1)) {
275                          sprintf(errmsg,
276                          "bad column count near line %d in sensor file '%s'",
# Line 282 | Line 321 | init_ptable(
321          if (pvals != NULL)
322                  free((void *)pvals);
323          if (sfile == NULL || !*sfile) {
324 +                sensor = NULL;
325 +                sntp[0] = sntp[1] = 0;
326                  pvals = NULL;
327                  ntheta = nphi = 0;
328                  return;
# Line 299 | Line 340 | init_ptable(
340                  error(INTERNAL, errmsg);
341          }
342                                          /* compute boundary angles */
343 <        maxtheta = 1.5f*sensor[sntp[0]*(sntp[1]+1)] -
344 <                        0.5f*sensor[sntp[0]*sntp[1]];
343 >        maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
344 >        if (maxtheta > PI)
345 >                maxtheta = PI;
346          thdiv[0] = .0;
347          for (t = 1; t < sntp[0]; t++)
348 <                thdiv[t] = DEGREE/2.*(sensor[t*(sntp[1]+1)] +
349 <                                        sensor[(t+1)*(sntp[1]+1)]);
350 <        thdiv[sntp[0]] = maxtheta*DEGREE;
309 <        phdiv[0] = .0;
348 >                thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
349 >        thdiv[sntp[0]] = maxtheta;
350 >        phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
351          for (p = 1; p < sntp[1]; p++)
352 <                phdiv[p] = DEGREE/2.*(sensor[p] + sensor[p+1]);
353 <        phdiv[sntp[1]] = 2.*PI;
352 >                phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
353 >        phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
354                                          /* size our table */
355 <        tsize = 1. - cos(maxtheta*DEGREE);
355 >        tsize = 1. - cos(maxtheta);
356          psize = PI*tsize/maxtheta;
357          if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
358                  samptot = sntp[0]*sntp[1];
359 <        ntheta = (int)(sqrt(samptot*tsize/psize) + 0.5);
359 >        ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
360          if (ntheta > MAXNT)
361                  ntheta = MAXNT;
362          nphi = samptot/ntheta;
363 <        pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
363 >        pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
364          if (pvals == NULL)
365                  error(SYSTEM, "out of memory in init_ptable()");
366          gscale = .0;                    /* compute our inverse table */
367          for (i = 0; i < sntp[0]; i++) {
368 <                rowp = sensor + (i+1)*(sntp[1]+1) + 1;
369 <                rowsum[i] = 0.;
368 >                rowp = &s_val(i,0);
369 >                rowsum[i] = 1e-20;
370                  for (j = 0; j < sntp[1]; j++)
371                          rowsum[i] += *rowp++;
372                  rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
373                  rowomega[i] *= 2.*PI / (double)sntp[1];
374                  gscale += rowsum[i] * rowomega[i];
375          }
376 <        tvals[0] = .0f;
377 <        for (i = 1; i < ntheta; i++) {
376 >        if (gscale <= FTINY) {
377 >                sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
378 >                error(USER, errmsg);
379 >        }
380 >        for (i = 0; i < ntheta; i++) {
381                  prob = (double)i / (double)ntheta;
382                  for (t = 0; t < sntp[0]; t++)
383                          if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
# Line 343 | Line 387 | init_ptable(
387                  frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
388                  tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
389                                                  frac*cos(thdiv[t+1]) );
390 <                pvals[i*(nphi+1)] = .0f;
390 >                                /* offset b/c sensor values are centered */
391 >                if (t <= 0 || frac > 0.5)
392 >                        frac -= 0.5;
393 >                else if (t >= sntp[0]-1 || frac < 0.5) {
394 >                        frac += 0.5;
395 >                        --t;
396 >                }
397 >                pvals[i*(nphi+1)] = phdiv[0];
398                  for (j = 1; j < nphi; j++) {
399                          prob = (double)j / (double)nphi;
400 <                        rowp = sensor + t*(sntp[1]+1) + 1;
401 <                        rowp1 = rowp + sntp[1]+1;
402 <                        for (p = 0; p < sntp[1]; p++) {
403 <                                if ((prob -= (1.-frac)*rowp[p]/rowsum[t-1] +
404 <                                                frac*rowp1[p]/rowsum[t]) <= .0)
400 >                        rowp = &s_val(t,0);
401 >                        rowp1 = &s_val(t+1,0);
402 >                        for (p = 0; p < sntp[1]; p++)
403 >                                if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
404 >                                            frac*rowp1[p]/rowsum[t+1]) <= .0)
405                                          break;
406 <                                if (p >= sntp[1])
407 <                                        error(INTERNAL,
408 <                                            "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];
406 >                        if (p >= sntp[1]) {
407 >                                p = sntp[1] - 1;
408 >                                prob = .5;
409                          }
410 +                        frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
411 +                                        + frac*rowp1[p]/rowsum[t+1]);
412 +                        pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
413 +                                                frac1*phdiv[p+1];
414                  }
415 <                pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
415 >                pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
416          }
417 +        tvals[0] = .0f;
418          tvals[ntheta] = (float)tsize;
419   }
420  
# Line 418 | Line 470 | sens_val(
470          t = (int)(theta/maxtheta * sntp[0]);
471          p = (int)(phi*(1./360.) * sntp[1]);
472                          /* hack for non-uniform sensor grid */
473 <        while (t+1 < sntp[0] && theta >= sensor[(t+2)*(sntp[1]+1)])
473 >        while (t+1 < sntp[0] && theta >= s_theta(t+1))
474                  ++t;
475 <        while (t-1 >= 0 && theta < sensor[t*(sntp[1]+1)])
475 >        while (t-1 >= 0 && theta <= s_theta(t-1))
476                  --t;
477 <        while (p+1 < sntp[1] && phi >= sensor[p+2])
477 >        while (p+1 < sntp[1] && phi >= s_phi(p+1))
478                  ++p;
479 <        while (p-1 >= 0 && phi < sensor[p])
479 >        while (p-1 >= 0 && phi <= s_phi(p-1))
480                  --p;
481 <        return(sensor[t*(sntp[1]+1) + p + 1]);
481 >        return(s_val(t,p));
482   }
483  
484 + /* Print origin and direction */
485 + static void
486 + print_ray(
487 +        FVECT rorg,
488 +        FVECT rdir
489 + )
490 + {
491 +        printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
492 +                        rorg[0], rorg[1], rorg[2],
493 +                        rdir[0], rdir[1], rdir[2]);
494 + }
495 +
496   /* Compute sensor output */
497   static void
498   comp_sensor(
# Line 441 | Line 505 | comp_sensor(
505          int     nt, np;
506          COLOR   vsum;
507          RAY     rr;
508 +        double  sf;
509          int     i, j;
510                                                  /* set view */
511          ourview.type = VT_ANG;
# Line 451 | Line 516 | comp_sensor(
516                  error(USER, err);
517                                                  /* assign probability table */
518          init_ptable(sfile);
519 <                                                /* do Monte Carlo sampling */
519 >                                                /* stratified MC sampling */
520          setcolor(vsum, .0f, .0f, .0f);
521          nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
522          np = nsamps/nt;
523 <        VCOPY(rr.rorg, ourview.vp);
459 <        rr.rmax = .0;
523 >        sf = gscale/nsamps;
524          for (i = 0; i < nt; i++)
525 <                for (j =0; j < np; j++) {
526 <                        get_direc(rr.rdir, (i+frandom())/nt,
527 <                                        (j + frandom())/np);
525 >                for (j = 0; j < np; j++) {
526 >                        VCOPY(rr.rorg, ourview.vp);
527 >                        get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
528 >                        if (ourview.vfore > FTINY)
529 >                                VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
530 >                        if (!ray_pnprocs) {
531 >                                print_ray(rr.rorg, rr.rdir);
532 >                                continue;
533 >                        }
534 >                        rr.rmax = .0;
535                          rayorigin(&rr, PRIMARY, NULL, NULL);
536 +                        scalecolor(rr.rcoef, sf);
537                          if (ray_pqueue(&rr) == 1)
538                                  addcolor(vsum, rr.rcol);
539                  }
540 <                                                /* finish MC calculation */
541 <        while (ray_presult(&rr, 0) > 0)
542 <                addcolor(vsum, rr.rcol);
543 <        scalecolor(vsum, gscale/(nt*np));
544 <                                                /* compute direct component */
540 >                                                /* remaining rays pure MC */
541 >        for (i = nsamps - nt*np; i-- > 0; ) {
542 >                VCOPY(rr.rorg, ourview.vp);
543 >                get_direc(rr.rdir, frandom(), frandom());
544 >                if (ourview.vfore > FTINY)
545 >                        VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
546 >                if (!ray_pnprocs) {
547 >                        print_ray(rr.rorg, rr.rdir);
548 >                        continue;
549 >                }
550 >                rr.rmax = .0;
551 >                rayorigin(&rr, PRIMARY, NULL, NULL);
552 >                scalecolor(rr.rcoef, sf);
553 >                if (ray_pqueue(&rr) == 1)
554 >                        addcolor(vsum, rr.rcol);
555 >        }
556 >        if (!ray_pnprocs)                       /* just printing rays */
557 >                return;
558 >                                                /* scale partial result */
559 >        scalecolor(vsum, sf);
560 >                                                /* add direct component */
561          for (i = ndirs; i-- > 0; ) {
562                  SRCINDEX        si;
563                  initsrcindex(&si);
564                  while (srcray(&rr, NULL, &si)) {
565 <                        double  d = sens_val(rr.rdir);
566 <                        if (d <= FTINY)
565 >                        sf = sens_val(rr.rdir);
566 >                        if (sf <= FTINY)
567                                  continue;
568 <                        d *= si.dom/ndirs;
569 <                        scalecolor(rr.rcoef, d);
568 >                        sf *= si.dom/ndirs;
569 >                        scalecolor(rr.rcoef, sf);
570                          if (ray_pqueue(&rr) == 1) {
571                                  multcolor(rr.rcol, rr.rcoef);
572                                  addcolor(vsum, rr.rcol);
573                          }
574                  }
575          }
576 <                                                /* finish direct calculation */
576 >                                                /* finish our calculation */
577          while (ray_presult(&rr, 0) > 0) {
578                  multcolor(rr.rcol, rr.rcoef);
579                  addcolor(vsum, rr.rcol);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines