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.6 by greg, Sat Dec 13 00:44:05 2008 UTC vs.
Revision 2.14 by greg, Tue Sep 30 23:14:34 2014 UTC

# Line 34 | Line 34 | int            nprocs = 1;     /* number of rendering processes
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 71 | Line 71 | print_defaults()               /* print out default parameters */
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 109 | Line 120 | main(
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]);   /* process a sensor file */
127                          continue;
# Line 183 | Line 195 | main(
195                  }
196                  i += rval;
197          }
198 <        if (!ray_pnprocs)
198 >        if (sensor == NULL)
199                  error(USER, i<argc ? "missing sensor file" : "missing octree");
200          quit(0);
201   }
# Line 195 | Line 207 | load_sensor(
207          char    *sfile
208   )
209   {
210 +        int     warnedneg;
211          char    linebuf[8192];
212          int     nelem = 1000;
213          float   *sarr = (float *)malloc(sizeof(float)*nelem);
# Line 223 | Line 236 | load_sensor(
236                  cp = fskip(cp);
237                  if (cp == NULL)
238                          break;
239 +                if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]) {
240 +                        sprintf(errmsg,
241 +                "Phi values not monotinically increasing in sensor file '%s'",
242 +                                        sfile);
243 +                        error(USER, errmsg);
244 +                }
245                  ++ntp[1];
246          }
247 +        warnedneg = 0;
248          ntp[0] = 0;                             /* get thetas + data */
249          while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
250                  ++ntp[0];
# Line 242 | Line 262 | load_sensor(
262                          cp = fskip(cp);
263                          if (cp == NULL)
264                                  break;
265 +                        if (i && sarr[i] < .0) {
266 +                                if (!warnedneg++) {
267 +                                        sprintf(errmsg,
268 +                "Negative value(s) in sensor file '%s' (ignored)\n", sfile);
269 +                                        error(WARNING, errmsg);
270 +                                }
271 +                                sarr[i] = .0;
272 +                        }
273                          ++i;
274                  }
275                  if (i == ntp[0]*(ntp[1]+1))
276                          break;
277 +                if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
278 +                                        sarr[(ntp[0]-1)*(ntp[1]+1)]) {
279 +                        sprintf(errmsg,
280 +                "Theta values not monotinically increasing in sensor file '%s'",
281 +                                        sfile);
282 +                        error(USER, errmsg);
283 +                }
284                  if (i != (ntp[0]+1)*(ntp[1]+1)) {
285                          sprintf(errmsg,
286                          "bad column count near line %d in sensor file '%s'",
# Line 315 | Line 350 | init_ptable(
350                  error(INTERNAL, errmsg);
351          }
352                                          /* compute boundary angles */
353 <        maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
353 >        maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
354 >        if (maxtheta > PI)
355 >                maxtheta = PI;
356          thdiv[0] = .0;
357          for (t = 1; t < sntp[0]; t++)
358                  thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
359 <        thdiv[sntp[0]] = maxtheta*DEGREE;
360 <        phdiv[0] = .0;
359 >        thdiv[sntp[0]] = maxtheta;
360 >        phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
361          for (p = 1; p < sntp[1]; p++)
362                  phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
363 <        phdiv[sntp[1]] = 2.*PI;
363 >        phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
364                                          /* size our table */
365 <        tsize = 1. - cos(maxtheta*DEGREE);
366 <        psize = PI*tsize/(maxtheta*DEGREE);
365 >        tsize = 1. - cos(maxtheta);
366 >        psize = PI*tsize/maxtheta;
367          if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
368                  samptot = sntp[0]*sntp[1];
369          ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
370          if (ntheta > MAXNT)
371                  ntheta = MAXNT;
372          nphi = samptot/ntheta;
373 <        pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
373 >        pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
374          if (pvals == NULL)
375                  error(SYSTEM, "out of memory in init_ptable()");
376          gscale = .0;                    /* compute our inverse table */
377          for (i = 0; i < sntp[0]; i++) {
378                  rowp = &s_val(i,0);
379 <                rowsum[i] = 0.;
379 >                rowsum[i] = 1e-20;
380                  for (j = 0; j < sntp[1]; j++)
381                          rowsum[i] += *rowp++;
382                  rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
383                  rowomega[i] *= 2.*PI / (double)sntp[1];
384                  gscale += rowsum[i] * rowomega[i];
385          }
386 +        if (gscale <= FTINY) {
387 +                sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
388 +                error(USER, errmsg);
389 +        }
390          for (i = 0; i < ntheta; i++) {
391                  prob = (double)i / (double)ntheta;
392                  for (t = 0; t < sntp[0]; t++)
# Line 357 | Line 398 | init_ptable(
398                  tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
399                                                  frac*cos(thdiv[t+1]) );
400                                  /* offset b/c sensor values are centered */
401 <                if (t <= 0 || frac > 0.5)
401 >                if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) {
402                          frac -= 0.5;
403 <                else if (t >= sntp[0]-1 || frac < 0.5) {
403 >                } else {
404                          frac += 0.5;
405                          --t;
406                  }
407 <                pvals[i*(nphi+1)] = .0f;
407 >                pvals[i*(nphi+1)] = phdiv[0];
408                  for (j = 1; j < nphi; j++) {
409                          prob = (double)j / (double)nphi;
410                          rowp = &s_val(t,0);
411                          rowp1 = &s_val(t+1,0);
412 <                        for (p = 0; p < sntp[1]; p++) {
412 >                        for (p = 0; p < sntp[1]; p++)
413                                  if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
414                                              frac*rowp1[p]/rowsum[t+1]) <= .0)
415                                          break;
416 <                                if (p >= sntp[1])
417 <                                        error(INTERNAL,
418 <                                            "code error 2 in init_ptable()");
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];
416 >                        if (p >= sntp[1]) {     /* should never happen? */
417 >                                p = sntp[1] - 1;
418 >                                prob = .5;
419                          }
420 +                        frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
421 +                                        + frac*rowp1[p]/rowsum[t+1]);
422 +                        pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
423 +                                                frac1*phdiv[p+1];
424                  }
425 <                pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
425 >                pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
426          }
427          tvals[0] = .0f;
428          tvals[ntheta] = (float)tsize;
# Line 429 | Line 470 | sens_val(
470          int     t, p;
471          
472          dv[2] = DOT(dvec, ourview.vdir);
473 <        theta = (float)((1./DEGREE) * acos(dv[2]));
473 >        theta = acos(dv[2]);
474          if (theta >= maxtheta)
475                  return(.0f);
476          dv[0] = DOT(dvec, ourview.hvec);
477          dv[1] = DOT(dvec, ourview.vvec);
478 <        phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
479 <        while (phi < .0f) phi += 360.f;
478 >        phi = atan2(-dv[0], dv[1]);
479 >        while (phi < .0f) phi += (float)(2.*PI);
480          t = (int)(theta/maxtheta * sntp[0]);
481 <        p = (int)(phi*(1./360.) * sntp[1]);
481 >        p = (int)(phi*(1./(2.*PI)) * sntp[1]);
482                          /* hack for non-uniform sensor grid */
483 +        theta *= (float)(1./DEGREE);
484 +        phi *= (float)(1./DEGREE);
485          while (t+1 < sntp[0] && theta >= s_theta(t+1))
486                  ++t;
487          while (t-1 >= 0 && theta <= s_theta(t-1))
# Line 450 | Line 493 | sens_val(
493          return(s_val(t,p));
494   }
495  
496 + /* Print origin and direction */
497 + static void
498 + print_ray(
499 +        FVECT rorg,
500 +        FVECT rdir
501 + )
502 + {
503 +        printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
504 +                        rorg[0], rorg[1], rorg[2],
505 +                        rdir[0], rdir[1], rdir[2]);
506 + }
507 +
508   /* Compute sensor output */
509   static void
510   comp_sensor(
# Line 484 | Line 539 | comp_sensor(
539                          get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
540                          if (ourview.vfore > FTINY)
541                                  VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
542 +                        if (!ray_pnprocs) {
543 +                                print_ray(rr.rorg, rr.rdir);
544 +                                continue;
545 +                        }
546                          rr.rmax = .0;
547 <                        rayorigin(&rr, PRIMARY, NULL, NULL);
547 >                        rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
548                          scalecolor(rr.rcoef, sf);
549                          if (ray_pqueue(&rr) == 1)
550                                  addcolor(vsum, rr.rcol);
# Line 496 | Line 555 | comp_sensor(
555                  get_direc(rr.rdir, frandom(), frandom());
556                  if (ourview.vfore > FTINY)
557                          VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
558 +                if (!ray_pnprocs) {
559 +                        print_ray(rr.rorg, rr.rdir);
560 +                        continue;
561 +                }
562                  rr.rmax = .0;
563 <                rayorigin(&rr, PRIMARY, NULL, NULL);
563 >                rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
564                  scalecolor(rr.rcoef, sf);
565                  if (ray_pqueue(&rr) == 1)
566                          addcolor(vsum, rr.rcol);
567          }
568 +        if (!ray_pnprocs)                       /* just printing rays */
569 +                return;
570                                                  /* scale partial result */
571          scalecolor(vsum, sf);
572                                                  /* add direct component */

Diff Legend

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