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.9 by greg, Tue May 17 19:34:36 2011 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 235 | Line 235 | load_sensor(
235                  cp = fskip(cp);
236                  if (cp == NULL)
237                          break;
238 +                if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]) {
239 +                        sprintf(errmsg,
240 +                "Phi values not monotinically increasing in sensor file '%s'",
241 +                                        sfile);
242 +                        error(USER, errmsg);
243 +                }
244                  ++ntp[1];
245          }
246          ntp[0] = 0;                             /* get thetas + data */
# Line 258 | Line 264 | load_sensor(
264                  }
265                  if (i == ntp[0]*(ntp[1]+1))
266                          break;
267 +                if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
268 +                                        sarr[(ntp[0]-1)*(ntp[1]+1)]) {
269 +                        sprintf(errmsg,
270 +                "Theta values not monotinically increasing in sensor file '%s'",
271 +                                        sfile);
272 +                        error(USER, errmsg);
273 +                }
274                  if (i != (ntp[0]+1)*(ntp[1]+1)) {
275                          sprintf(errmsg,
276                          "bad column count near line %d in sensor file '%s'",
# Line 327 | Line 340 | init_ptable(
340                  error(INTERNAL, errmsg);
341          }
342                                          /* compute boundary angles */
343 <        maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
343 >        maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
344 >        if (maxtheta > PI)
345 >                maxtheta = PI;
346          thdiv[0] = .0;
347          for (t = 1; t < sntp[0]; t++)
348                  thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
349 <        thdiv[sntp[0]] = maxtheta*DEGREE;
350 <        phdiv[0] = .0;
349 >        thdiv[sntp[0]] = maxtheta;
350 >        phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
351          for (p = 1; p < sntp[1]; p++)
352                  phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
353 <        phdiv[sntp[1]] = 2.*PI;
353 >        phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
354                                          /* size our table */
355 <        tsize = 1. - cos(maxtheta*DEGREE);
356 <        psize = PI*tsize/(maxtheta*DEGREE);
355 >        tsize = 1. - cos(maxtheta);
356 >        psize = PI*tsize/maxtheta;
357          if (sntp[0]*sntp[1] < samptot)  /* don't overdo resolution */
358                  samptot = sntp[0]*sntp[1];
359          ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
360          if (ntheta > MAXNT)
361                  ntheta = MAXNT;
362          nphi = samptot/ntheta;
363 <        pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
363 >        pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
364          if (pvals == NULL)
365                  error(SYSTEM, "out of memory in init_ptable()");
366          gscale = .0;                    /* compute our inverse table */
367          for (i = 0; i < sntp[0]; i++) {
368                  rowp = &s_val(i,0);
369 <                rowsum[i] = 0.;
369 >                rowsum[i] = 1e-20;
370                  for (j = 0; j < sntp[1]; j++)
371                          rowsum[i] += *rowp++;
372                  rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
373                  rowomega[i] *= 2.*PI / (double)sntp[1];
374                  gscale += rowsum[i] * rowomega[i];
375          }
376 +        if (gscale <= FTINY) {
377 +                sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
378 +                error(USER, errmsg);
379 +        }
380          for (i = 0; i < ntheta; i++) {
381                  prob = (double)i / (double)ntheta;
382                  for (t = 0; t < sntp[0]; t++)
# Line 375 | Line 394 | init_ptable(
394                          frac += 0.5;
395                          --t;
396                  }
397 <                pvals[i*(nphi+1)] = .0f;
397 >                pvals[i*(nphi+1)] = phdiv[0];
398                  for (j = 1; j < nphi; j++) {
399                          prob = (double)j / (double)nphi;
400                          rowp = &s_val(t,0);
401                          rowp1 = &s_val(t+1,0);
402 <                        for (p = 0; p < sntp[1]; p++) {
402 >                        for (p = 0; p < sntp[1]; p++)
403                                  if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
404                                              frac*rowp1[p]/rowsum[t+1]) <= .0)
405                                          break;
406 <                                if (p >= sntp[1])
407 <                                        error(INTERNAL,
408 <                                            "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];
406 >                        if (p >= sntp[1]) {
407 >                                p = sntp[1] - 1;
408 >                                prob = .5;
409                          }
410 +                        frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
411 +                                        + frac*rowp1[p]/rowsum[t+1]);
412 +                        pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
413 +                                                frac1*phdiv[p+1];
414                  }
415 <                pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
415 >                pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
416          }
417          tvals[0] = .0f;
418          tvals[ntheta] = (float)tsize;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines