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.8 by greg, Sun Sep 26 15:41:46 2010 UTC vs.
Revision 2.21 by greg, Tue Mar 15 18:05:03 2022 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] */
# Line 26 | Line 27 | VIEW           ourview =  {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,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 */
30 < 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 59 | 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);
62 –        /* 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 77 | Line 76 | quit(ec)                       /* make sure exit is called */
76   int     ec;
77   {
78          if (ray_pnprocs > 0)    /* close children if any */
79 <                ray_pclose(0);          
79 >                ray_pclose(0);
80 >        else if (ray_pnprocs < 0)
81 >                _exit(ec);      /* avoid flush in child */
82          exit(ec);
83   }
84  
# Line 129 | Line 130 | main(
130                  if (argv[i][1] == 'r') {        /* sampling options */
131                          if (argv[i][2] == 'd')
132                                  nsamps = atol(argv[++i]);
132 –                        else if (argv[i][2] == 's')
133 –                                nssamps = atol(argv[++i]);
133                          else {
134                                  sprintf(errmsg, "bad option at '%s'", argv[i]);
135                                  error(USER, errmsg);
# Line 207 | Line 206 | load_sensor(
206          char    *sfile
207   )
208   {
209 +        int     warnedneg;
210          char    linebuf[8192];
211 +        int     last_pos_val = 0;
212          int     nelem = 1000;
213          float   *sarr = (float *)malloc(sizeof(float)*nelem);
214          FILE    *fp;
# Line 235 | 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]]+FTINY) {
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 254 | Line 262 | load_sensor(
262                          cp = fskip(cp);
263                          if (cp == NULL)
264                                  break;
265 +                        if (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 +                        } else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1))
273 +                                last_pos_val = i;
274                          ++i;
275                  }
276 <                if (i == ntp[0]*(ntp[1]+1))
276 >                if (i == ntp[0]*(ntp[1]+1))     /* empty line? */
277                          break;
278 +                if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
279 +                                        sarr[(ntp[0]-1)*(ntp[1]+1)]) {
280 +                        sprintf(errmsg,
281 +                "Theta values not monotinically increasing in sensor file '%s'",
282 +                                        sfile);
283 +                        error(USER, errmsg);
284 +                }
285                  if (i != (ntp[0]+1)*(ntp[1]+1)) {
286                          sprintf(errmsg,
287                          "bad column count near line %d in sensor file '%s'",
# Line 265 | Line 289 | load_sensor(
289                          error(USER, errmsg);
290                  }
291          }
292 <        nelem = i;
292 >                                                /* truncate zero region */
293 >        ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1;
294 >        nelem = (ntp[0]+1)*(ntp[1]+1);
295          fclose(fp);
296          errmsg[0] = '\0';                       /* sanity checks */
297 <        if (ntp[0] <= 0)
298 <                sprintf(errmsg, "no data in sensor file '%s'", sfile);
297 >        if (!last_pos_val)
298 >                sprintf(errmsg, "no positive sensor values in file '%s'", sfile);
299          else if (fabs(sarr[ntp[1]+1]) > FTINY)
300                  sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
301                                  sfile);
302          else if (fabs(sarr[1]) > FTINY)
303                  sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
304                                  sfile);
305 <        else if (sarr[ntp[1]] <= FTINY)
305 >        else if (sarr[ntp[1]] < 270.-FTINY)
306                  sprintf(errmsg,
307 <                        "maximum phi must be positive in sensor file '%s'",
307 >                        "maximum phi must be 270 or greater in sensor file '%s'",
308                                  sfile);
309 <        else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
309 >        else if (sarr[ntp[1]] >= 360.-FTINY)
310                  sprintf(errmsg,
311 <                        "maximum theta must be positive in sensor file '%s'",
311 >                        "maximum phi must be less than 360 in sensor file '%s'",
312                                  sfile);
313          if (errmsg[0])
314                  error(USER, errmsg);
# Line 295 | Line 321 | init_ptable(
321          char    *sfile
322   )
323   {
324 <        int     samptot = nsamps;
324 >        long    samptot = nsamps;
325          float   *rowp, *rowp1;
326          double  rowsum[MAXNT], rowomega[MAXNT];
327          double  thdiv[MAXNT+1], phdiv[MAXNP+1];
# Line 327 | Line 353 | init_ptable(
353                  error(INTERNAL, errmsg);
354          }
355                                          /* compute boundary angles */
356 <        maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
356 >        maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
357 >        if (maxtheta > PI)
358 >                maxtheta = PI;
359          thdiv[0] = .0;
360          for (t = 1; t < sntp[0]; t++)
361                  thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
362 <        thdiv[sntp[0]] = maxtheta*DEGREE;
363 <        phdiv[0] = .0;
362 >        thdiv[sntp[0]] = maxtheta;
363 >        phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
364          for (p = 1; p < sntp[1]; p++)
365                  phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
366 <        phdiv[sntp[1]] = 2.*PI;
366 >        phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
367                                          /* size our table */
368 <        tsize = 1. - cos(maxtheta*DEGREE);
369 <        psize = PI*tsize/(maxtheta*DEGREE);
368 >        tsize = 1. - cos(maxtheta);
369 >        psize = PI*tsize/maxtheta;
370          if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
371                  samptot = sntp[0]*sntp[1];
372 <        ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
372 >        ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1;
373          if (ntheta > MAXNT)
374                  ntheta = MAXNT;
375          nphi = samptot/ntheta;
376 <        pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
376 >        pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
377          if (pvals == NULL)
378                  error(SYSTEM, "out of memory in init_ptable()");
379          gscale = .0;                    /* compute our inverse table */
380          for (i = 0; i < sntp[0]; i++) {
381                  rowp = &s_val(i,0);
382 <                rowsum[i] = 0.;
382 >                rowsum[i] = 1e-20;
383                  for (j = 0; j < sntp[1]; j++)
384                          rowsum[i] += *rowp++;
385                  rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
386                  rowomega[i] *= 2.*PI / (double)sntp[1];
387                  gscale += rowsum[i] * rowomega[i];
388          }
389 +        if (gscale <= FTINY) {
390 +                sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
391 +                error(USER, errmsg);
392 +        }
393          for (i = 0; i < ntheta; i++) {
394                  prob = (double)i / (double)ntheta;
395                  for (t = 0; t < sntp[0]; t++)
# Line 369 | Line 401 | init_ptable(
401                  tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
402                                                  frac*cos(thdiv[t+1]) );
403                                  /* offset b/c sensor values are centered */
404 <                if (t <= 0 || frac > 0.5)
404 >                if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) {
405                          frac -= 0.5;
406 <                else if (t >= sntp[0]-1 || frac < 0.5) {
406 >                } else {
407                          frac += 0.5;
408                          --t;
409                  }
410 <                pvals[i*(nphi+1)] = .0f;
410 >                pvals[i*(nphi+1)] = phdiv[0];
411                  for (j = 1; j < nphi; j++) {
412                          prob = (double)j / (double)nphi;
413                          rowp = &s_val(t,0);
414                          rowp1 = &s_val(t+1,0);
415 <                        for (p = 0; p < sntp[1]; p++) {
415 >                        for (p = 0; p < sntp[1]; p++)
416                                  if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
417                                              frac*rowp1[p]/rowsum[t+1]) <= .0)
418                                          break;
419 <                                if (p >= sntp[1])
420 <                                        error(INTERNAL,
421 <                                            "code error 2 in init_ptable()");
390 <                                frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
391 <                                                + frac*rowp1[p]/rowsum[t+1]);
392 <                                pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
393 <                                                        frac1*phdiv[p+1];
419 >                        if (p >= sntp[1]) {     /* should never happen? */
420 >                                p = sntp[1] - 1;
421 >                                prob = .5;
422                          }
423 +                        frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
424 +                                        + frac*rowp1[p]/rowsum[t+1]);
425 +                        pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
426 +                                                frac1*phdiv[p+1];
427                  }
428 <                pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
428 >                pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
429          }
430 +                                                /* duplicate final row */
431 +        memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1),
432 +                                sizeof(*pvals)*(nphi+1));
433          tvals[0] = .0f;
434          tvals[ntheta] = (float)tsize;
435   }
# Line 441 | Line 476 | sens_val(
476          int     t, p;
477          
478          dv[2] = DOT(dvec, ourview.vdir);
479 <        theta = (float)((1./DEGREE) * acos(dv[2]));
479 >        theta = acos(dv[2]);
480          if (theta >= maxtheta)
481                  return(.0f);
482          dv[0] = DOT(dvec, ourview.hvec);
483          dv[1] = DOT(dvec, ourview.vvec);
484 <        phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
485 <        while (phi < .0f) phi += 360.f;
484 >        phi = atan2(-dv[0], dv[1]);
485 >        while (phi < .0f) phi += (float)(2.*PI);
486          t = (int)(theta/maxtheta * sntp[0]);
487 <        p = (int)(phi*(1./360.) * sntp[1]);
487 >        p = (int)(phi*(1./(2.*PI)) * sntp[1]);
488                          /* hack for non-uniform sensor grid */
489 +        theta *= (float)(1./DEGREE);
490 +        phi *= (float)(1./DEGREE);
491          while (t+1 < sntp[0] && theta >= s_theta(t+1))
492                  ++t;
493          while (t-1 >= 0 && theta <= s_theta(t-1))
# Line 513 | Line 550 | comp_sensor(
550                                  continue;
551                          }
552                          rr.rmax = .0;
553 <                        rayorigin(&rr, PRIMARY, NULL, NULL);
553 >                        rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
554                          scalecolor(rr.rcoef, sf);
555                          if (ray_pqueue(&rr) == 1)
556                                  addcolor(vsum, rr.rcol);
# Line 529 | Line 566 | comp_sensor(
566                          continue;
567                  }
568                  rr.rmax = .0;
569 <                rayorigin(&rr, PRIMARY, NULL, NULL);
569 >                rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL);
570                  scalecolor(rr.rcoef, sf);
571                  if (ray_pqueue(&rr) == 1)
572                          addcolor(vsum, rr.rcol);

Diff Legend

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