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

Comparing src/util/rsensor.c (file contents):
Revision 2.4 by greg, Fri Feb 22 21:52:10 2008 UTC vs.
Revision 2.25 by greg, Tue Jun 3 21:31:51 2025 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 "func.h"
15   #include "view.h"
16   #include "random.h"
17  
18   #define DEGREE          (PI/180.)
19  
20 < #define MAXNT           180     /* maximum number of theta divisions */
20 > #define MAXNT           181     /* maximum number of theta divisions */
21   #define MAXNP           360     /* maximum number of phi divisions */
22  
21 – 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 */
28 < unsigned long   nssamps = 9000; /* number of super-samples */
30 > long            nsamps = 10000; /* desired number of initial 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 */
# Line 57 | Line 59 | print_defaults()               /* print out default parameters */
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);
60 –        /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
62          printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
63          printf("-vp %f %f %f\t# view point\n",
64                          ourview.vp[0], ourview.vp[1], ourview.vp[2]);
# Line 69 | Line 70 | print_defaults()               /* print out default parameters */
70          print_rdefaults();
71   }
72  
73 +
74 + void
75 + quit(ec)                        /* make sure exit is called */
76 + int     ec;
77 + {
78 +        if (ray_pnprocs > 0)    /* close children if any */
79 +                ray_pclose(0);
80 +        else if (ray_pnprocs < 0)
81 +                _exit(ec);      /* avoid flush in child */
82 +        exit(ec);
83 + }
84 +
85 +
86   int
87   main(
88          int     argc,
# Line 78 | Line 92 | main(
92          int     doheader = 1;
93          int     optwarn = 0;
94          int     i, rval;
95 <
96 <        progname = argv[0];
95 >                                /* set global progname */
96 >        fixargv0(argv[0]);
97                                  /* set up rendering defaults */
98          rand_samp = 1;
99 <        dstrsrc = 0.5;
99 >        dstrsrc = 0.65;
100          srcsizerat = 0.1;
101          directrelay = 3;
102          ambounce = 1;
103          maxdepth = -10;
104 +                                /* initialize calcomp routines */
105 +        initfunc();
106                                  /* get options from command line */
107          for (i = 1; i < argc; i++) {
108                  while ((rval = expandarg(&argc, &argv, i)) > 0)
# Line 101 | Line 117 | main(
117                          if (!ray_pnprocs) {
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]);   /* process a sensor file */
130                          continue;
# Line 114 | Line 132 | main(
132                  if (argv[i][1] == 'r') {        /* sampling options */
133                          if (argv[i][2] == 'd')
134                                  nsamps = atol(argv[++i]);
117 –                        else if (argv[i][2] == 's')
118 –                                nssamps = atol(argv[++i]);
135                          else {
136                                  sprintf(errmsg, "bad option at '%s'", argv[i]);
137                                  error(USER, errmsg);
# Line 180 | Line 196 | main(
196                  }
197                  i += rval;
198          }
199 <        if (!ray_pnprocs)
199 >        if (sensor == NULL)
200                  error(USER, i<argc ? "missing sensor file" : "missing octree");
201          quit(0);
202   }
# Line 192 | Line 208 | load_sensor(
208          char    *sfile
209   )
210   {
211 +        int     warnedneg;
212          char    linebuf[8192];
213 +        int     last_pos_val = 0;
214          int     nelem = 1000;
215          float   *sarr = (float *)malloc(sizeof(float)*nelem);
216          FILE    *fp;
# Line 220 | Line 238 | load_sensor(
238                  cp = fskip(cp);
239                  if (cp == NULL)
240                          break;
241 +                if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]+FTINY) {
242 +                        sprintf(errmsg,
243 +                "Phi values not monotinically increasing in sensor file '%s'",
244 +                                        sfile);
245 +                        error(USER, errmsg);
246 +                }
247                  ++ntp[1];
248          }
249 +        warnedneg = 0;
250          ntp[0] = 0;                             /* get thetas + data */
251          while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
252                  ++ntp[0];
# Line 239 | Line 264 | load_sensor(
264                          cp = fskip(cp);
265                          if (cp == NULL)
266                                  break;
267 +                        if (sarr[i] < .0) {
268 +                                if (!warnedneg++) {
269 +                                        sprintf(errmsg,
270 +                "Negative value(s) in sensor file '%s' (ignored)\n", sfile);
271 +                                        error(WARNING, errmsg);
272 +                                }
273 +                                sarr[i] = .0;
274 +                        } else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1))
275 +                                last_pos_val = i;
276                          ++i;
277                  }
278 <                if (i == ntp[0]*(ntp[1]+1))
278 >                if (i == ntp[0]*(ntp[1]+1))     /* empty line? */
279                          break;
280 +                if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
281 +                                        sarr[(ntp[0]-1)*(ntp[1]+1)]) {
282 +                        sprintf(errmsg,
283 +                "Theta values not monotinically increasing in sensor file '%s'",
284 +                                        sfile);
285 +                        error(USER, errmsg);
286 +                }
287                  if (i != (ntp[0]+1)*(ntp[1]+1)) {
288                          sprintf(errmsg,
289                          "bad column count near line %d in sensor file '%s'",
# Line 250 | Line 291 | load_sensor(
291                          error(USER, errmsg);
292                  }
293          }
294 <        nelem = i;
294 >                                                /* truncate zero region */
295 >        ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1;
296 >        nelem = (ntp[0]+1)*(ntp[1]+1);
297          fclose(fp);
298          errmsg[0] = '\0';                       /* sanity checks */
299 <        if (ntp[0] <= 0)
300 <                sprintf(errmsg, "no data in sensor file '%s'", sfile);
299 >        if (!last_pos_val)
300 >                sprintf(errmsg, "no positive sensor values in file '%s'", sfile);
301          else if (fabs(sarr[ntp[1]+1]) > FTINY)
302                  sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
303                                  sfile);
304          else if (fabs(sarr[1]) > FTINY)
305                  sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
306                                  sfile);
307 <        else if (sarr[ntp[1]] <= FTINY)
307 >        else if (sarr[ntp[1]] < 270.-FTINY)
308                  sprintf(errmsg,
309 <                        "maximum phi must be positive in sensor file '%s'",
309 >                        "maximum phi must be 270 or greater in sensor file '%s'",
310                                  sfile);
311 <        else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
311 >        else if (sarr[ntp[1]] >= 360.-FTINY)
312                  sprintf(errmsg,
313 <                        "maximum theta must be positive in sensor file '%s'",
313 >                        "maximum phi must be less than 360 in sensor file '%s'",
314                                  sfile);
315          if (errmsg[0])
316                  error(USER, errmsg);
# Line 280 | Line 323 | init_ptable(
323          char    *sfile
324   )
325   {
326 <        int     samptot = nsamps;
326 >        long    samptot = nsamps;
327          float   *rowp, *rowp1;
328          double  rowsum[MAXNT], rowomega[MAXNT];
329          double  thdiv[MAXNT+1], phdiv[MAXNP+1];
# Line 312 | Line 355 | init_ptable(
355                  error(INTERNAL, errmsg);
356          }
357                                          /* compute boundary angles */
358 <        maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
358 >        maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
359 >        if (maxtheta > PI)
360 >                maxtheta = PI;
361          thdiv[0] = .0;
362          for (t = 1; t < sntp[0]; t++)
363                  thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
364 <        thdiv[sntp[0]] = maxtheta*DEGREE;
365 <        phdiv[0] = .0;
364 >        thdiv[sntp[0]] = maxtheta;
365 >        phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
366          for (p = 1; p < sntp[1]; p++)
367                  phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
368 <        phdiv[sntp[1]] = 2.*PI;
368 >        phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
369                                          /* size our table */
370 <        tsize = 1. - cos(maxtheta*DEGREE);
371 <        psize = PI*tsize/(maxtheta*DEGREE);
370 >        tsize = 1. - cos(maxtheta);
371 >        psize = PI*tsize/maxtheta;
372          if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
373                  samptot = sntp[0]*sntp[1];
374 <        ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
374 >        ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1;
375          if (ntheta > MAXNT)
376                  ntheta = MAXNT;
377          nphi = samptot/ntheta;
378 <        pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
378 >        pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
379          if (pvals == NULL)
380                  error(SYSTEM, "out of memory in init_ptable()");
381          gscale = .0;                    /* compute our inverse table */
382          for (i = 0; i < sntp[0]; i++) {
383                  rowp = &s_val(i,0);
384 <                rowsum[i] = 0.;
384 >                rowsum[i] = 1e-20;
385                  for (j = 0; j < sntp[1]; j++)
386                          rowsum[i] += *rowp++;
387                  rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
388                  rowomega[i] *= 2.*PI / (double)sntp[1];
389                  gscale += rowsum[i] * rowomega[i];
390          }
391 +        if (gscale <= FTINY) {
392 +                sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
393 +                error(USER, errmsg);
394 +        }
395          for (i = 0; i < ntheta; i++) {
396                  prob = (double)i / (double)ntheta;
397                  for (t = 0; t < sntp[0]; t++)
# Line 354 | Line 403 | init_ptable(
403                  tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
404                                                  frac*cos(thdiv[t+1]) );
405                                  /* offset b/c sensor values are centered */
406 <                if (t <= 0 || frac > 0.5)
406 >                if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) {
407                          frac -= 0.5;
408 <                else if (t >= sntp[0]-1 || frac < 0.5) {
408 >                } else {
409                          frac += 0.5;
410                          --t;
411                  }
412 <                pvals[i*(nphi+1)] = .0f;
412 >                pvals[i*(nphi+1)] = phdiv[0];
413                  for (j = 1; j < nphi; j++) {
414                          prob = (double)j / (double)nphi;
415                          rowp = &s_val(t,0);
416                          rowp1 = &s_val(t+1,0);
417 <                        for (p = 0; p < sntp[1]; p++) {
417 >                        for (p = 0; p < sntp[1]; p++)
418                                  if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
419                                              frac*rowp1[p]/rowsum[t+1]) <= .0)
420                                          break;
421 <                                if (p >= sntp[1])
422 <                                        error(INTERNAL,
423 <                                            "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];
421 >                        if (p >= sntp[1]) {     /* should never happen? */
422 >                                p = sntp[1] - 1;
423 >                                prob = .5;
424                          }
425 +                        frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
426 +                                        + frac*rowp1[p]/rowsum[t+1]);
427 +                        pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
428 +                                                frac1*phdiv[p+1];
429                  }
430 <                pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
430 >                pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
431          }
432 +                                                /* duplicate final row */
433 +        memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1),
434 +                                sizeof(*pvals)*(nphi+1));
435          tvals[0] = .0f;
436          tvals[ntheta] = (float)tsize;
437   }
# Line 426 | Line 478 | sens_val(
478          int     t, p;
479          
480          dv[2] = DOT(dvec, ourview.vdir);
481 <        theta = (float)((1./DEGREE) * acos(dv[2]));
481 >        theta = acos(dv[2]);
482          if (theta >= maxtheta)
483                  return(.0f);
484          dv[0] = DOT(dvec, ourview.hvec);
485          dv[1] = DOT(dvec, ourview.vvec);
486 <        phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
487 <        while (phi < .0f) phi += 360.f;
486 >        phi = atan2(-dv[0], dv[1]);
487 >        while (phi < .0f) phi += (float)(2.*PI);
488          t = (int)(theta/maxtheta * sntp[0]);
489 <        p = (int)(phi*(1./360.) * sntp[1]);
489 >        p = (int)(phi*(1./(2.*PI)) * sntp[1]);
490                          /* hack for non-uniform sensor grid */
491 +        theta *= (float)(1./DEGREE);
492 +        phi *= (float)(1./DEGREE);
493          while (t+1 < sntp[0] && theta >= s_theta(t+1))
494                  ++t;
495          while (t-1 >= 0 && theta <= s_theta(t-1))
# Line 447 | Line 501 | sens_val(
501          return(s_val(t,p));
502   }
503  
504 + /* Print origin and direction */
505 + static void
506 + print_ray(
507 +        FVECT rorg,
508 +        FVECT rdir
509 + )
510 + {
511 +        printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
512 +                        rorg[0], rorg[1], rorg[2],
513 +                        rdir[0], rdir[1], rdir[2]);
514 + }
515 +
516   /* Compute sensor output */
517   static void
518   comp_sensor(
# Line 457 | Line 523 | comp_sensor(
523                                  ndsamps > 0 ? 1 : 0;
524          char    *err;
525          int     nt, np;
526 <        COLOR   vsum;
526 >        SCOLOR  vsum;
527          RAY     rr;
528 +        double  sf;
529          int     i, j;
530                                                  /* set view */
531          ourview.type = VT_ANG;
# Line 469 | Line 536 | comp_sensor(
536                  error(USER, err);
537                                                  /* assign probability table */
538          init_ptable(sfile);
539 <                                                /* do Monte Carlo sampling */
540 <        setcolor(vsum, .0f, .0f, .0f);
539 >                                                /* stratified MC sampling */
540 >        scolorblack(vsum);
541          nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
542          np = nsamps/nt;
543 <        VCOPY(rr.rorg, ourview.vp);
477 <        rr.rmax = .0;
543 >        sf = gscale/nsamps;
544          for (i = 0; i < nt; i++)
545                  for (j = 0; j < np; j++) {
546 +                        VCOPY(rr.rorg, ourview.vp);
547                          get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
548 <                        rayorigin(&rr, PRIMARY, NULL, NULL);
548 >                        if (ourview.vfore > FTINY)
549 >                                VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
550 >                        if (!ray_pnprocs) {
551 >                                print_ray(rr.rorg, rr.rdir);
552 >                                continue;
553 >                        }
554 >                        rr.rmax = .0;
555 >                        rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
556 >                        scalescolor(rr.rcoef, sf);
557                          if (ray_pqueue(&rr) == 1)
558 <                                addcolor(vsum, rr.rcol);
558 >                                saddscolor(vsum, rr.rcol);
559                  }
560 <                                                /* finish MC calculation */
561 <        while (ray_presult(&rr, 0) > 0)
562 <                addcolor(vsum, rr.rcol);
563 <        scalecolor(vsum, gscale/(nt*np));
564 <                                                /* compute direct component */
560 >                                                /* remaining rays pure MC */
561 >        for (i = nsamps - nt*np; i-- > 0; ) {
562 >                VCOPY(rr.rorg, ourview.vp);
563 >                get_direc(rr.rdir, frandom(), frandom());
564 >                if (ourview.vfore > FTINY)
565 >                        VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
566 >                if (!ray_pnprocs) {
567 >                        print_ray(rr.rorg, rr.rdir);
568 >                        continue;
569 >                }
570 >                rr.rmax = .0;
571 >                rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
572 >                scalescolor(rr.rcoef, sf);
573 >                if (ray_pqueue(&rr) == 1)
574 >                        saddscolor(vsum, rr.rcol);
575 >        }
576 >        if (!ray_pnprocs)                       /* just printing rays */
577 >                return;
578 >                                                /* scale partial result */
579 >        scalescolor(vsum, sf);
580 >                                                /* add direct component */
581          for (i = ndirs; i-- > 0; ) {
582                  SRCINDEX        si;
583                  initsrcindex(&si);
584                  while (srcray(&rr, NULL, &si)) {
585 <                        double  d = sens_val(rr.rdir);
586 <                        if (d <= FTINY)
585 >                        sf = sens_val(rr.rdir);
586 >                        if (sf <= FTINY)
587                                  continue;
588 <                        d *= si.dom/ndirs;
589 <                        scalecolor(rr.rcoef, d);
588 >                        sf *= si.dom/ndirs;
589 >                        scalescolor(rr.rcoef, sf);
590                          if (ray_pqueue(&rr) == 1) {
591 <                                multcolor(rr.rcol, rr.rcoef);
592 <                                addcolor(vsum, rr.rcol);
591 >                                smultscolor(rr.rcol, rr.rcoef);
592 >                                saddscolor(vsum, rr.rcol);
593                          }
594                  }
595          }
596 <                                                /* finish direct calculation */
597 <        while (ray_presult(&rr, 0) > 0) {
598 <                multcolor(rr.rcol, rr.rcoef);
508 <                addcolor(vsum, rr.rcol);
596 >        while (ray_presult(&rr, 0) > 0) {       /* finish our calculation */
597 >                smultscolor(rr.rcol, rr.rcoef);
598 >                saddscolor(vsum, rr.rcol);
599          }
600 <                                                /* print our result */
601 <        printf("%.4e %.4e %.4e\n", colval(vsum,RED),
602 <                                colval(vsum,GRN), colval(vsum,BLU));
600 >        for (i = 0; i < NCSAMP; i++)            /* print our result */
601 >                printf(" %.4e", vsum[i]);
602 >        fputc('\n', stdout);
603   }

Diff Legend

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