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.2 by greg, Thu Feb 21 20:18:25 2008 UTC vs.
Revision 2.20 by greg, Fri Jul 24 17:09:33 2020 UTC

# Line 9 | Line 9 | static const char RCSid[] = "$Id$";
9   */
10  
11   #include "ray.h"
12 + #include "platform.h"
13   #include "source.h"
14   #include "view.h"
15   #include "random.h"
16  
17   #define DEGREE          (PI/180.)
18  
19 < #define MAXNT           180     /* maximum number of theta divisions */
19 > #define MAXNT           181     /* maximum number of theta divisions */
20   #define MAXNP           360     /* maximum number of phi divisions */
21  
22   extern char     *progname;      /* global argv[0] */
23   extern int      nowarn;         /* don't report warnings? */
24  
25                                  /* current sensor's perspective */
26 < VIEW            ourview = STDVIEW;
26 > VIEW            ourview =  {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,0.,0.},
27 >                                1.,180.,180.,0.,0.,0.,0.,
28 >                                {0.,0.,0.},{0.,0.,0.},0.,0.};
29  
30   unsigned long   nsamps = 10000; /* desired number of initial samples */
31   unsigned long   nssamps = 9000; /* number of super-samples */
# Line 32 | Line 35 | int            nprocs = 1;     /* number of rendering processes
35   float           *sensor = NULL; /* current sensor data */
36   int             sntp[2];        /* number of sensor theta and phi angles */
37   float           maxtheta;       /* maximum theta value for this sensor */
38 < float           tvals[MAXNT+1]; /* theta values (1-D table of 1-cos(t)) */
39 < float           *pvals = NULL;  /* phi values (2-D table in radians) */
38 > float           tvals[MAXNT+1]; /* theta prob. values (1-D table of 1-cos(t)) */
39 > float           *pvals = NULL;  /* phi prob. values (2-D table in radians) */
40   int             ntheta = 0;     /* polar angle divisions */
41   int             nphi = 0;       /* azimuthal angle divisions */
42   double          gscale = 1.;    /* global scaling value */
# Line 45 | Line 48 | double         gscale = 1.;    /* global scaling value */
48   static void     comp_sensor(char *sfile);
49  
50   static void
51 < print_defaults()
51 > over_options()                  /* overriding options */
52   {
53 +        directvis = (ndsamps <= 0);
54 +        do_irrad = 0;
55 + }
56 +
57 + static void
58 + print_defaults()                /* print out default parameters */
59 + {
60 +        over_options();
61          printf("-n %-9d\t\t\t# number of processes\n", nprocs);
62          printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
63          /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
# Line 61 | Line 72 | print_defaults()
72          print_rdefaults();
73   }
74  
75 +
76 + void
77 + quit(ec)                        /* make sure exit is called */
78 + int     ec;
79 + {
80 +        if (ray_pnprocs > 0)    /* close children if any */
81 +                ray_pclose(0);
82 +        else if (ray_pnprocs < 0)
83 +                _exit(ec);      /* avoid flush in child */
84 +        exit(ec);
85 + }
86 +
87 +
88   int
89   main(
90          int     argc,
# Line 68 | Line 92 | main(
92   )
93   {
94          int     doheader = 1;
95 +        int     optwarn = 0;
96          int     i, rval;
97  
98          progname = argv[0];
99                                  /* set up rendering defaults */
100          rand_samp = 1;
101 <        dstrsrc = 0.5;
101 >        dstrsrc = 0.65;
102          srcsizerat = 0.1;
103          directrelay = 3;
104          ambounce = 1;
105          maxdepth = -10;
81 –                                /* just asking defaults? */
82 –        if (argc == 2 && !strcmp(argv[1], "-defaults")) {
83 –                print_defaults();
84 –                return(0);
85 –        }
86 –                                /* check octree */
87 –        if (argc < 2 || argv[argc-1][0] == '-')
88 –                error(USER, "missing octree argument");
106                                  /* get options from command line */
107 <        for (i = 1; i < argc-1; i++) {
107 >        for (i = 1; i < argc; i++) {
108                  while ((rval = expandarg(&argc, &argv, i)) > 0)
109                          ;
110                  if (rval < 0) {
111                          sprintf(errmsg, "cannot expand '%s'", argv[i]);
112                          error(SYSTEM, errmsg);
113                  }
114 <                if (argv[i][0] != '-') {        /* process a sensor file */
114 >                if (argv[i][0] != '-') {
115 >                        if (i >= argc-1)
116 >                                break;          /* final octree argument */
117                          if (!ray_pnprocs) {
118 <                                                /* overriding options */
100 <                                directvis = (ndsamps <= 0);
101 <                                do_irrad = 0;
118 >                                over_options();
119                                  if (doheader) { /* print header */
120 +                                        newheader("RADIANCE", stdout);
121                                          printargs(argc, argv, stdout);
122                                          fputformat("ascii", stdout);
123                                          putchar('\n');
124                                  }
125                                                  /* start process(es) */
126 <                                ray_pinit(argv[argc-1], nprocs);
126 >                                if (strcmp(argv[argc-1], "."))
127 >                                        ray_pinit(argv[argc-1], nprocs);
128                          }
129 <                        comp_sensor(argv[i]);
129 >                        comp_sensor(argv[i]);   /* process a sensor file */
130                          continue;
131                  }
132                  if (argv[i][1] == 'r') {        /* sampling options */
# Line 150 | Line 169 | main(
169                          sprintf(errmsg, "bad view option at '%s'", argv[i]);
170                          error(USER, errmsg);
171                  }
172 <                if (!strcmp(argv[i], "-w")) {   /* turn off warnings */
173 <                        nowarn = 1;
172 >                if (!strcmp(argv[i], "-w")) {   /* toggle warnings */
173 >                        nowarn = !nowarn;
174                          continue;
175                  }
176                  if (ray_pnprocs) {
177 <                        error(WARNING,
177 >                        if (!optwarn++)
178 >                                error(WARNING,
179                          "rendering options should appear before first sensor");
180                  } else if (!strcmp(argv[i], "-defaults")) {
181                          print_defaults();
# Line 178 | Line 198 | main(
198                  }
199                  i += rval;
200          }
201 +        if (sensor == NULL)
202 +                error(USER, i<argc ? "missing sensor file" : "missing octree");
203          quit(0);
204   }
205  
# Line 188 | Line 210 | load_sensor(
210          char    *sfile
211   )
212   {
213 +        int     warnedneg;
214          char    linebuf[8192];
215 +        int     last_pos_val = 0;
216          int     nelem = 1000;
217          float   *sarr = (float *)malloc(sizeof(float)*nelem);
218          FILE    *fp;
# Line 216 | Line 240 | load_sensor(
240                  cp = fskip(cp);
241                  if (cp == NULL)
242                          break;
243 +                if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]+FTINY) {
244 +                        sprintf(errmsg,
245 +                "Phi values not monotinically increasing in sensor file '%s'",
246 +                                        sfile);
247 +                        error(USER, errmsg);
248 +                }
249                  ++ntp[1];
250          }
251 +        warnedneg = 0;
252          ntp[0] = 0;                             /* get thetas + data */
253          while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
254                  ++ntp[0];
# Line 235 | Line 266 | load_sensor(
266                          cp = fskip(cp);
267                          if (cp == NULL)
268                                  break;
269 +                        if (sarr[i] < .0) {
270 +                                if (!warnedneg++) {
271 +                                        sprintf(errmsg,
272 +                "Negative value(s) in sensor file '%s' (ignored)\n", sfile);
273 +                                        error(WARNING, errmsg);
274 +                                }
275 +                                sarr[i] = .0;
276 +                        } else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1))
277 +                                last_pos_val = i;
278                          ++i;
279                  }
280 <                if (i == ntp[0]*(ntp[1]+1))
280 >                if (i == ntp[0]*(ntp[1]+1))     /* empty line? */
281                          break;
282 +                if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
283 +                                        sarr[(ntp[0]-1)*(ntp[1]+1)]) {
284 +                        sprintf(errmsg,
285 +                "Theta values not monotinically increasing in sensor file '%s'",
286 +                                        sfile);
287 +                        error(USER, errmsg);
288 +                }
289                  if (i != (ntp[0]+1)*(ntp[1]+1)) {
290                          sprintf(errmsg,
291                          "bad column count near line %d in sensor file '%s'",
# Line 246 | Line 293 | load_sensor(
293                          error(USER, errmsg);
294                  }
295          }
296 <        nelem = i;
296 >                                                /* truncate zero region */
297 >        ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1;
298 >        nelem = (ntp[0]+1)*(ntp[1]+1);
299          fclose(fp);
300          errmsg[0] = '\0';                       /* sanity checks */
301 <        if (ntp[0] <= 0)
302 <                sprintf(errmsg, "no data in sensor file '%s'", sfile);
301 >        if (!last_pos_val)
302 >                sprintf(errmsg, "no positive sensor values in file '%s'", sfile);
303          else if (fabs(sarr[ntp[1]+1]) > FTINY)
304                  sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
305                                  sfile);
306          else if (fabs(sarr[1]) > FTINY)
307                  sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
308                                  sfile);
309 <        else if (sarr[ntp[1]] <= FTINY)
309 >        else if (sarr[ntp[1]] < 270.-FTINY)
310                  sprintf(errmsg,
311 <                        "maximum phi must be positive in sensor file '%s'",
311 >                        "maximum phi must be 270 or greater in sensor file '%s'",
312                                  sfile);
313 <        else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
313 >        else if (sarr[ntp[1]] >= 360.-FTINY)
314                  sprintf(errmsg,
315 <                        "maximum theta must be positive in sensor file '%s'",
315 >                        "maximum phi must be less than 360 in sensor file '%s'",
316                                  sfile);
317          if (errmsg[0])
318                  error(USER, errmsg);
# Line 308 | Line 357 | init_ptable(
357                  error(INTERNAL, errmsg);
358          }
359                                          /* compute boundary angles */
360 <        maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
360 >        maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
361 >        if (maxtheta > PI)
362 >                maxtheta = PI;
363          thdiv[0] = .0;
364          for (t = 1; t < sntp[0]; t++)
365                  thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
366 <        thdiv[sntp[0]] = maxtheta*DEGREE;
367 <        phdiv[0] = .0;
366 >        thdiv[sntp[0]] = maxtheta;
367 >        phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
368          for (p = 1; p < sntp[1]; p++)
369                  phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
370 <        phdiv[sntp[1]] = 2.*PI;
370 >        phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
371                                          /* size our table */
372 <        tsize = 1. - cos(maxtheta*DEGREE);
373 <        psize = PI*tsize/(maxtheta*DEGREE);
372 >        tsize = 1. - cos(maxtheta);
373 >        psize = PI*tsize/maxtheta;
374          if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
375                  samptot = sntp[0]*sntp[1];
376 <        ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
376 >        ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1;
377          if (ntheta > MAXNT)
378                  ntheta = MAXNT;
379          nphi = samptot/ntheta;
380 <        pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
380 >        pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
381          if (pvals == NULL)
382                  error(SYSTEM, "out of memory in init_ptable()");
383          gscale = .0;                    /* compute our inverse table */
384          for (i = 0; i < sntp[0]; i++) {
385                  rowp = &s_val(i,0);
386 <                rowsum[i] = 0.;
386 >                rowsum[i] = 1e-20;
387                  for (j = 0; j < sntp[1]; j++)
388                          rowsum[i] += *rowp++;
389                  rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
390                  rowomega[i] *= 2.*PI / (double)sntp[1];
391                  gscale += rowsum[i] * rowomega[i];
392          }
393 +        if (gscale <= FTINY) {
394 +                sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
395 +                error(USER, errmsg);
396 +        }
397          for (i = 0; i < ntheta; i++) {
398                  prob = (double)i / (double)ntheta;
399                  for (t = 0; t < sntp[0]; t++)
# Line 350 | Line 405 | init_ptable(
405                  tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
406                                                  frac*cos(thdiv[t+1]) );
407                                  /* offset b/c sensor values are centered */
408 <                if (t <= 0 || frac > 0.5)
408 >                if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) {
409                          frac -= 0.5;
410 <                else if (t >= sntp[0]-1 || frac < 0.5) {
410 >                } else {
411                          frac += 0.5;
412                          --t;
413                  }
414 <                pvals[i*(nphi+1)] = .0f;
414 >                pvals[i*(nphi+1)] = phdiv[0];
415                  for (j = 1; j < nphi; j++) {
416                          prob = (double)j / (double)nphi;
417                          rowp = &s_val(t,0);
418                          rowp1 = &s_val(t+1,0);
419 <                        for (p = 0; p < sntp[1]; p++) {
419 >                        for (p = 0; p < sntp[1]; p++)
420                                  if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
421                                              frac*rowp1[p]/rowsum[t+1]) <= .0)
422                                          break;
423 <                                if (p >= sntp[1])
424 <                                        error(INTERNAL,
425 <                                            "code error 2 in init_ptable()");
371 <                                frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
372 <                                                + frac*rowp1[p]/rowsum[t+1]);
373 <                                if (p <= 0 || frac1 > 0.5)
374 <                                        frac1 -= 0.5;
375 <                                else if (p >= sntp[1]-1 || frac1 < 0.5) {
376 <                                        frac1 += 0.5;
377 <                                        --p;
378 <                                }
379 <                                pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
380 <                                                        frac1*phdiv[p+1];
423 >                        if (p >= sntp[1]) {     /* should never happen? */
424 >                                p = sntp[1] - 1;
425 >                                prob = .5;
426                          }
427 +                        frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
428 +                                        + frac*rowp1[p]/rowsum[t+1]);
429 +                        pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
430 +                                                frac1*phdiv[p+1];
431                  }
432 <                pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
432 >                pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
433          }
434 +                                                /* duplicate final row */
435 +        memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1),
436 +                                sizeof(*pvals)*(nphi+1));
437          tvals[0] = .0f;
438          tvals[ntheta] = (float)tsize;
439   }
# Line 428 | Line 480 | sens_val(
480          int     t, p;
481          
482          dv[2] = DOT(dvec, ourview.vdir);
483 <        theta = (float)((1./DEGREE) * acos(dv[2]));
483 >        theta = acos(dv[2]);
484          if (theta >= maxtheta)
485                  return(.0f);
486          dv[0] = DOT(dvec, ourview.hvec);
487          dv[1] = DOT(dvec, ourview.vvec);
488 <        phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
489 <        while (phi < .0f) phi += 360.f;
488 >        phi = atan2(-dv[0], dv[1]);
489 >        while (phi < .0f) phi += (float)(2.*PI);
490          t = (int)(theta/maxtheta * sntp[0]);
491 <        p = (int)(phi*(1./360.) * sntp[1]);
491 >        p = (int)(phi*(1./(2.*PI)) * sntp[1]);
492                          /* hack for non-uniform sensor grid */
493 +        theta *= (float)(1./DEGREE);
494 +        phi *= (float)(1./DEGREE);
495          while (t+1 < sntp[0] && theta >= s_theta(t+1))
496                  ++t;
497          while (t-1 >= 0 && theta <= s_theta(t-1))
# Line 449 | Line 503 | sens_val(
503          return(s_val(t,p));
504   }
505  
506 + /* Print origin and direction */
507 + static void
508 + print_ray(
509 +        FVECT rorg,
510 +        FVECT rdir
511 + )
512 + {
513 +        printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
514 +                        rorg[0], rorg[1], rorg[2],
515 +                        rdir[0], rdir[1], rdir[2]);
516 + }
517 +
518   /* Compute sensor output */
519   static void
520   comp_sensor(
# Line 461 | Line 527 | comp_sensor(
527          int     nt, np;
528          COLOR   vsum;
529          RAY     rr;
530 +        double  sf;
531          int     i, j;
532                                                  /* set view */
533          ourview.type = VT_ANG;
# Line 471 | Line 538 | comp_sensor(
538                  error(USER, err);
539                                                  /* assign probability table */
540          init_ptable(sfile);
541 <                                                /* do Monte Carlo sampling */
541 >                                                /* stratified MC sampling */
542          setcolor(vsum, .0f, .0f, .0f);
543          nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
544          np = nsamps/nt;
545 <        VCOPY(rr.rorg, ourview.vp);
479 <        rr.rmax = .0;
545 >        sf = gscale/nsamps;
546          for (i = 0; i < nt; i++)
547                  for (j = 0; j < np; j++) {
548 +                        VCOPY(rr.rorg, ourview.vp);
549                          get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
550 <                        rayorigin(&rr, PRIMARY, NULL, NULL);
550 >                        if (ourview.vfore > FTINY)
551 >                                VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
552 >                        if (!ray_pnprocs) {
553 >                                print_ray(rr.rorg, rr.rdir);
554 >                                continue;
555 >                        }
556 >                        rr.rmax = .0;
557 >                        rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
558 >                        scalecolor(rr.rcoef, sf);
559                          if (ray_pqueue(&rr) == 1)
560                                  addcolor(vsum, rr.rcol);
561                  }
562 <                                                /* finish MC calculation */
563 <        while (ray_presult(&rr, 0) > 0)
564 <                addcolor(vsum, rr.rcol);
565 <        scalecolor(vsum, gscale/(nt*np));
566 <                                                /* compute direct component */
562 >                                                /* remaining rays pure MC */
563 >        for (i = nsamps - nt*np; i-- > 0; ) {
564 >                VCOPY(rr.rorg, ourview.vp);
565 >                get_direc(rr.rdir, frandom(), frandom());
566 >                if (ourview.vfore > FTINY)
567 >                        VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
568 >                if (!ray_pnprocs) {
569 >                        print_ray(rr.rorg, rr.rdir);
570 >                        continue;
571 >                }
572 >                rr.rmax = .0;
573 >                rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
574 >                scalecolor(rr.rcoef, sf);
575 >                if (ray_pqueue(&rr) == 1)
576 >                        addcolor(vsum, rr.rcol);
577 >        }
578 >        if (!ray_pnprocs)                       /* just printing rays */
579 >                return;
580 >                                                /* scale partial result */
581 >        scalecolor(vsum, sf);
582 >                                                /* add direct component */
583          for (i = ndirs; i-- > 0; ) {
584                  SRCINDEX        si;
585                  initsrcindex(&si);
586                  while (srcray(&rr, NULL, &si)) {
587 <                        double  d = sens_val(rr.rdir);
588 <                        if (d <= FTINY)
587 >                        sf = sens_val(rr.rdir);
588 >                        if (sf <= FTINY)
589                                  continue;
590 <                        d *= si.dom/ndirs;
591 <                        scalecolor(rr.rcoef, d);
590 >                        sf *= si.dom/ndirs;
591 >                        scalecolor(rr.rcoef, sf);
592                          if (ray_pqueue(&rr) == 1) {
593                                  multcolor(rr.rcol, rr.rcoef);
594                                  addcolor(vsum, rr.rcol);
595                          }
596                  }
597          }
598 <                                                /* finish direct calculation */
598 >                                                /* finish our calculation */
599          while (ray_presult(&rr, 0) > 0) {
600                  multcolor(rr.rcol, rr.rcoef);
601                  addcolor(vsum, rr.rcol);

Diff Legend

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