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.6 by greg, Sat Dec 13 00:44:05 2008 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 */
# Line 38 | Line 40 | 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 64 | Line 78 | main(
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 <        dstrsrc = 0.25;
86 >        rand_samp = 1;
87 >        dstrsrc = 0.65;
88 >        srcsizerat = 0.1;
89          directrelay = 3;
90          ambounce = 1;
91 <                                /* 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");
91 >        maxdepth = -10;
92                                  /* get options from command line */
93 <        for (i = 1; i < argc-1; i++) {
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] != '-') {        /* process a sensor file */
100 >                if (argv[i][0] != '-') {
101 >                        if (i >= argc-1)
102 >                                break;          /* final octree argument */
103                          if (!ray_pnprocs) {
104 <                                                /* overriding options */
93 <                                directvis = (ndsamps <= 0);
94 <                                do_irrad = 0;
104 >                                over_options();
105                                  if (doheader) { /* print header */
106 +                                        newheader("RADIANCE", stdout);
107                                          printargs(argc, argv, stdout);
108                                          fputformat("ascii", stdout);
109                                          putchar('\n');
# Line 100 | Line 111 | main(
111                                                  /* start process(es) */
112                                  ray_pinit(argv[argc-1], nprocs);
113                          }
114 <                        comp_sensor(argv[i]);
114 >                        comp_sensor(argv[i]);   /* process a sensor file */
115                          continue;
116                  }
117                  if (argv[i][1] == 'r') {        /* sampling options */
# Line 143 | Line 154 | main(
154                          sprintf(errmsg, "bad view option at '%s'", argv[i]);
155                          error(USER, errmsg);
156                  }
157 <                if (!strcmp(argv[i], "-w")) {   /* turn off warnings */
158 <                        nowarn = 1;
157 >                if (!strcmp(argv[i], "-w")) {   /* toggle warnings */
158 >                        nowarn = !nowarn;
159                          continue;
160                  }
161                  if (ray_pnprocs) {
162 <                        error(WARNING,
162 >                        if (!optwarn++)
163 >                                error(WARNING,
164                          "rendering options should appear before first sensor");
165                  } else if (!strcmp(argv[i], "-defaults")) {
166                          print_defaults();
# Line 171 | Line 183 | main(
183                  }
184                  i += rval;
185          }
186 +        if (!ray_pnprocs)
187 +                error(USER, i<argc ? "missing sensor file" : "missing octree");
188          quit(0);
189   }
190  
# Line 282 | Line 296 | init_ptable(
296          if (pvals != NULL)
297                  free((void *)pvals);
298          if (sfile == NULL || !*sfile) {
299 +                sensor = NULL;
300 +                sntp[0] = sntp[1] = 0;
301                  pvals = NULL;
302                  ntheta = nphi = 0;
303                  return;
# Line 299 | Line 315 | init_ptable(
315                  error(INTERNAL, errmsg);
316          }
317                                          /* compute boundary angles */
318 <        maxtheta = 1.5f*sensor[sntp[0]*(sntp[1]+1)] -
303 <                        0.5f*sensor[sntp[0]*sntp[1]];
318 >        maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
319          thdiv[0] = .0;
320          for (t = 1; t < sntp[0]; t++)
321 <                thdiv[t] = DEGREE/2.*(sensor[t*(sntp[1]+1)] +
307 <                                        sensor[(t+1)*(sntp[1]+1)]);
321 >                thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
322          thdiv[sntp[0]] = maxtheta*DEGREE;
323          phdiv[0] = .0;
324          for (p = 1; p < sntp[1]; p++)
325 <                phdiv[p] = DEGREE/2.*(sensor[p] + sensor[p+1]);
325 >                phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
326          phdiv[sntp[1]] = 2.*PI;
327                                          /* size our table */
328          tsize = 1. - cos(maxtheta*DEGREE);
329 <        psize = PI*tsize/maxtheta;
329 >        psize = PI*tsize/(maxtheta*DEGREE);
330          if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
331                  samptot = sntp[0]*sntp[1];
332 <        ntheta = (int)(sqrt(samptot*tsize/psize) + 0.5);
332 >        ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
333          if (ntheta > MAXNT)
334                  ntheta = MAXNT;
335          nphi = samptot/ntheta;
# Line 324 | Line 338 | init_ptable(
338                  error(SYSTEM, "out of memory in init_ptable()");
339          gscale = .0;                    /* compute our inverse table */
340          for (i = 0; i < sntp[0]; i++) {
341 <                rowp = sensor + (i+1)*(sntp[1]+1) + 1;
341 >                rowp = &s_val(i,0);
342                  rowsum[i] = 0.;
343                  for (j = 0; j < sntp[1]; j++)
344                          rowsum[i] += *rowp++;
# Line 332 | Line 346 | init_ptable(
346                  rowomega[i] *= 2.*PI / (double)sntp[1];
347                  gscale += rowsum[i] * rowomega[i];
348          }
349 <        tvals[0] = .0f;
336 <        for (i = 1; i < ntheta; i++) {
349 >        for (i = 0; i < ntheta; i++) {
350                  prob = (double)i / (double)ntheta;
351                  for (t = 0; t < sntp[0]; t++)
352                          if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
# Line 343 | Line 356 | init_ptable(
356                  frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
357                  tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
358                                                  frac*cos(thdiv[t+1]) );
359 +                                /* offset b/c sensor values are centered */
360 +                if (t <= 0 || frac > 0.5)
361 +                        frac -= 0.5;
362 +                else if (t >= sntp[0]-1 || frac < 0.5) {
363 +                        frac += 0.5;
364 +                        --t;
365 +                }
366                  pvals[i*(nphi+1)] = .0f;
367                  for (j = 1; j < nphi; j++) {
368                          prob = (double)j / (double)nphi;
369 <                        rowp = sensor + t*(sntp[1]+1) + 1;
370 <                        rowp1 = rowp + sntp[1]+1;
369 >                        rowp = &s_val(t,0);
370 >                        rowp1 = &s_val(t+1,0);
371                          for (p = 0; p < sntp[1]; p++) {
372 <                                if ((prob -= (1.-frac)*rowp[p]/rowsum[t-1] +
373 <                                                frac*rowp1[p]/rowsum[t]) <= .0)
372 >                                if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
373 >                                            frac*rowp1[p]/rowsum[t+1]) <= .0)
374                                          break;
375                                  if (p >= sntp[1])
376                                          error(INTERNAL,
377                                              "code error 2 in init_ptable()");
378 <                                frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t-1]
379 <                                                + frac*rowp1[p]/rowsum[t]);
378 >                                frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
379 >                                                + frac*rowp1[p]/rowsum[t+1]);
380                                  pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
381                                                          frac1*phdiv[p+1];
382                          }
383                  }
384                  pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
385          }
386 +        tvals[0] = .0f;
387          tvals[ntheta] = (float)tsize;
388   }
389  
# Line 418 | Line 439 | sens_val(
439          t = (int)(theta/maxtheta * sntp[0]);
440          p = (int)(phi*(1./360.) * sntp[1]);
441                          /* hack for non-uniform sensor grid */
442 <        while (t+1 < sntp[0] && theta >= sensor[(t+2)*(sntp[1]+1)])
442 >        while (t+1 < sntp[0] && theta >= s_theta(t+1))
443                  ++t;
444 <        while (t-1 >= 0 && theta < sensor[t*(sntp[1]+1)])
444 >        while (t-1 >= 0 && theta <= s_theta(t-1))
445                  --t;
446 <        while (p+1 < sntp[1] && phi >= sensor[p+2])
446 >        while (p+1 < sntp[1] && phi >= s_phi(p+1))
447                  ++p;
448 <        while (p-1 >= 0 && phi < sensor[p])
448 >        while (p-1 >= 0 && phi <= s_phi(p-1))
449                  --p;
450 <        return(sensor[t*(sntp[1]+1) + p + 1]);
450 >        return(s_val(t,p));
451   }
452  
453   /* Compute sensor output */
# Line 441 | Line 462 | comp_sensor(
462          int     nt, np;
463          COLOR   vsum;
464          RAY     rr;
465 +        double  sf;
466          int     i, j;
467                                                  /* set view */
468          ourview.type = VT_ANG;
# Line 451 | Line 473 | comp_sensor(
473                  error(USER, err);
474                                                  /* assign probability table */
475          init_ptable(sfile);
476 <                                                /* do Monte Carlo sampling */
476 >                                                /* stratified MC sampling */
477          setcolor(vsum, .0f, .0f, .0f);
478          nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
479          np = nsamps/nt;
480 <        VCOPY(rr.rorg, ourview.vp);
459 <        rr.rmax = .0;
480 >        sf = gscale/nsamps;
481          for (i = 0; i < nt; i++)
482 <                for (j =0; j < np; j++) {
483 <                        get_direc(rr.rdir, (i+frandom())/nt,
484 <                                        (j + frandom())/np);
482 >                for (j = 0; j < np; j++) {
483 >                        VCOPY(rr.rorg, ourview.vp);
484 >                        get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
485 >                        if (ourview.vfore > FTINY)
486 >                                VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
487 >                        rr.rmax = .0;
488                          rayorigin(&rr, PRIMARY, NULL, NULL);
489 +                        scalecolor(rr.rcoef, sf);
490                          if (ray_pqueue(&rr) == 1)
491                                  addcolor(vsum, rr.rcol);
492                  }
493 <                                                /* finish MC calculation */
494 <        while (ray_presult(&rr, 0) > 0)
495 <                addcolor(vsum, rr.rcol);
496 <        scalecolor(vsum, gscale/(nt*np));
497 <                                                /* compute direct component */
493 >                                                /* remaining rays pure MC */
494 >        for (i = nsamps - nt*np; i-- > 0; ) {
495 >                VCOPY(rr.rorg, ourview.vp);
496 >                get_direc(rr.rdir, frandom(), frandom());
497 >                if (ourview.vfore > FTINY)
498 >                        VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
499 >                rr.rmax = .0;
500 >                rayorigin(&rr, PRIMARY, NULL, NULL);
501 >                scalecolor(rr.rcoef, sf);
502 >                if (ray_pqueue(&rr) == 1)
503 >                        addcolor(vsum, rr.rcol);
504 >        }
505 >                                                /* scale partial result */
506 >        scalecolor(vsum, sf);
507 >                                                /* add direct component */
508          for (i = ndirs; i-- > 0; ) {
509                  SRCINDEX        si;
510                  initsrcindex(&si);
511                  while (srcray(&rr, NULL, &si)) {
512 <                        double  d = sens_val(rr.rdir);
513 <                        if (d <= FTINY)
512 >                        sf = sens_val(rr.rdir);
513 >                        if (sf <= FTINY)
514                                  continue;
515 <                        d *= si.dom/ndirs;
516 <                        scalecolor(rr.rcoef, d);
515 >                        sf *= si.dom/ndirs;
516 >                        scalecolor(rr.rcoef, sf);
517                          if (ray_pqueue(&rr) == 1) {
518                                  multcolor(rr.rcol, rr.rcoef);
519                                  addcolor(vsum, rr.rcol);
520                          }
521                  }
522          }
523 <                                                /* finish direct calculation */
523 >                                                /* finish our calculation */
524          while (ray_presult(&rr, 0) > 0) {
525                  multcolor(rr.rcol, rr.rcoef);
526                  addcolor(vsum, rr.rcol);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines