Back to index

radiance  4R0+20100331
Defines | Functions | Variables
rsensor.c File Reference
#include "ray.h"
#include "source.h"
#include "view.h"
#include "random.h"

Go to the source code of this file.

Defines

#define DEGREE   (PI/180.)
#define MAXNT   180 /* maximum number of theta divisions */
#define MAXNP   360 /* maximum number of phi divisions */
#define s_theta(t)   sensor[(t+1)*(sntp[1]+1)]
#define s_phi(p)   sensor[(p)+1]
#define s_val(t, p)   sensor[(p)+1+(t+1)*(sntp[1]+1)]

Functions

static void comp_sensor (char *sfile)
static void over_options ()
static void print_defaults ()
void quit (int ec)
int main (int argc, char *argv[])
static float * load_sensor (int ntp[2], char *sfile)
static void init_ptable (char *sfile)
static void get_direc (FVECT dvec, double x, double y)
static float sens_val (FVECT dvec)

Variables

static const char RCSid [] = "$Id: rsensor.c,v 2.7 2009/12/12 19:01:00 greg Exp $"
char * progname
int nowarn
VIEW ourview
unsigned long nsamps = 10000
unsigned long nssamps = 9000
int ndsamps = 32
int nprocs = 1
float * sensor = NULL
int sntp [2]
float maxtheta
float tvals [MAXNT+1]
float * pvals = NULL
int ntheta = 0
int nphi = 0
double gscale = 1.

Define Documentation

#define DEGREE   (PI/180.)

Definition at line 16 of file rsensor.c.

#define MAXNP   360 /* maximum number of phi divisions */

Definition at line 19 of file rsensor.c.

#define MAXNT   180 /* maximum number of theta divisions */

Definition at line 18 of file rsensor.c.

#define s_phi (   p)    sensor[(p)+1]

Definition at line 44 of file rsensor.c.

#define s_theta (   t)    sensor[(t+1)*(sntp[1]+1)]

Definition at line 43 of file rsensor.c.

#define s_val (   t,
 
)    sensor[(p)+1+(t+1)*(sntp[1]+1)]

Definition at line 45 of file rsensor.c.


Function Documentation

static void comp_sensor ( char *  sfile) [static]

Definition at line 466 of file rsensor.c.

{
       int    ndirs = dstrsrc > FTINY ? ndsamps :
                            ndsamps > 0 ? 1 : 0;
       char   *err;
       int    nt, np;
       COLOR  vsum;
       RAY    rr;
       double sf;
       int    i, j;
                                          /* set view */
       ourview.type = VT_ANG;
       ourview.horiz = ourview.vert = 180.;
       ourview.hoff = ourview.voff = .0;
       err = setview(&ourview);
       if (err != NULL)
              error(USER, err);
                                          /* assign probability table */
       init_ptable(sfile);
                                          /* stratified MC sampling */
       setcolor(vsum, .0f, .0f, .0f);
       nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
       np = nsamps/nt;
       sf = gscale/nsamps;
       for (i = 0; i < nt; i++)
              for (j = 0; j < np; j++) {
                     VCOPY(rr.rorg, ourview.vp);
                     get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
                     if (ourview.vfore > FTINY)
                            VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
                     rr.rmax = .0;
                     rayorigin(&rr, PRIMARY, NULL, NULL);
                     scalecolor(rr.rcoef, sf);
                     if (ray_pqueue(&rr) == 1)
                            addcolor(vsum, rr.rcol);
              }
                                          /* remaining rays pure MC */
       for (i = nsamps - nt*np; i-- > 0; ) {
              VCOPY(rr.rorg, ourview.vp);
              get_direc(rr.rdir, frandom(), frandom());
              if (ourview.vfore > FTINY)
                     VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
              rr.rmax = .0;
              rayorigin(&rr, PRIMARY, NULL, NULL);
              scalecolor(rr.rcoef, sf);
              if (ray_pqueue(&rr) == 1)
                     addcolor(vsum, rr.rcol);
       }
                                          /* scale partial result */
       scalecolor(vsum, sf);
                                          /* add direct component */
       for (i = ndirs; i-- > 0; ) {
              SRCINDEX      si;
              initsrcindex(&si);
              while (srcray(&rr, NULL, &si)) {
                     sf = sens_val(rr.rdir);
                     if (sf <= FTINY)
                            continue;
                     sf *= si.dom/ndirs;
                     scalecolor(rr.rcoef, sf);
                     if (ray_pqueue(&rr) == 1) {
                            multcolor(rr.rcol, rr.rcoef);
                            addcolor(vsum, rr.rcol);
                     }
              }
       }
                                          /* finish our calculation */
       while (ray_presult(&rr, 0) > 0) {
              multcolor(rr.rcol, rr.rcoef);
              addcolor(vsum, rr.rcol);
       }
                                          /* print our result */
       printf("%.4e %.4e %.4e\n", colval(vsum,RED),
                            colval(vsum,GRN), colval(vsum,BLU));
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_direc ( FVECT  dvec,
double  x,
double  y 
) [static]

Definition at line 403 of file rsensor.c.

{
       double xfrac = x*ntheta;
       int    tndx = (int)xfrac;
       double yfrac = y*nphi;
       int    pndx = (int)yfrac;
       double rad, phi;
       FVECT  dv;
       int    i;

       xfrac -= (double)tndx;
       yfrac -= (double)pndx;
       pndx += tndx*(nphi+1);
       
       dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]);
       rad = sqrt(1. - dv[2]*dv[2]);
       phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1];
       dv[0] = -rad*sin(phi);
       dv[1] = rad*cos(phi);
       for (i = 3; i--; )
              dvec[i] = dv[0]*ourview.hvec[i] +
                            dv[1]*ourview.vvec[i] +
                            dv[2]*ourview.vdir[i] ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void init_ptable ( char *  sfile) [static]

Definition at line 293 of file rsensor.c.

{
       int    samptot = nsamps;
       float  *rowp, *rowp1;
       double rowsum[MAXNT], rowomega[MAXNT];
       double thdiv[MAXNT+1], phdiv[MAXNP+1];
       double tsize, psize;
       double prob, frac, frac1;
       int    i, j, t, p;
                                   /* free old table */
       if (sensor != NULL)
              free((void *)sensor);
       if (pvals != NULL)
              free((void *)pvals);
       if (sfile == NULL || !*sfile) {
              sensor = NULL;
              sntp[0] = sntp[1] = 0;
              pvals = NULL;
              ntheta = nphi = 0;
              return;
       }
                                   /* load sensor table */
       sensor = load_sensor(sntp, sfile);
       if (sntp[0] > MAXNT) {
              sprintf(errmsg, "Too many theta rows in sensor file '%s'",
                            sfile);
              error(INTERNAL, errmsg);
       }
       if (sntp[1] > MAXNP) {
              sprintf(errmsg, "Too many phi columns in sensor file '%s'",
                            sfile);
              error(INTERNAL, errmsg);
       }
                                   /* compute boundary angles */
       maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
       thdiv[0] = .0;
       for (t = 1; t < sntp[0]; t++)
              thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
       thdiv[sntp[0]] = maxtheta*DEGREE;
       phdiv[0] = .0;
       for (p = 1; p < sntp[1]; p++)
              phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
       phdiv[sntp[1]] = 2.*PI;
                                   /* size our table */
       tsize = 1. - cos(maxtheta*DEGREE);
       psize = PI*tsize/(maxtheta*DEGREE);
       if (sntp[0]*sntp[1] < samptot)     /* don't overdo resolution */
              samptot = sntp[0]*sntp[1];
       ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
       if (ntheta > MAXNT)
              ntheta = MAXNT;
       nphi = samptot/ntheta;
       pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
       if (pvals == NULL)
              error(SYSTEM, "out of memory in init_ptable()");
       gscale = .0;                /* compute our inverse table */
       for (i = 0; i < sntp[0]; i++) {
              rowp = &s_val(i,0);
              rowsum[i] = 0.;
              for (j = 0; j < sntp[1]; j++)
                     rowsum[i] += *rowp++;
              rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
              rowomega[i] *= 2.*PI / (double)sntp[1];
              gscale += rowsum[i] * rowomega[i];
       }
       for (i = 0; i < ntheta; i++) {
              prob = (double)i / (double)ntheta;
              for (t = 0; t < sntp[0]; t++)
                     if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
                            break;
              if (t >= sntp[0])
                     error(INTERNAL, "code error 1 in init_ptable()");
              frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
              tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
                                          frac*cos(thdiv[t+1]) );
                            /* offset b/c sensor values are centered */
              if (t <= 0 || frac > 0.5)
                     frac -= 0.5;
              else if (t >= sntp[0]-1 || frac < 0.5) {
                     frac += 0.5;
                     --t;
              }
              pvals[i*(nphi+1)] = .0f;
              for (j = 1; j < nphi; j++) {
                     prob = (double)j / (double)nphi;
                     rowp = &s_val(t,0);
                     rowp1 = &s_val(t+1,0);
                     for (p = 0; p < sntp[1]; p++) {
                            if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
                                       frac*rowp1[p]/rowsum[t+1]) <= .0)
                                   break;
                            if (p >= sntp[1])
                                   error(INTERNAL,
                                       "code error 2 in init_ptable()");
                            frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
                                          + frac*rowp1[p]/rowsum[t+1]);
                            pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
                                                 frac1*phdiv[p+1];
                     }
              }
              pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
       }
       tvals[0] = .0f;
       tvals[ntheta] = (float)tsize;
}

Here is the call graph for this function:

Here is the caller graph for this function:

static float* load_sensor ( int  ntp[2],
char *  sfile 
) [static]

Definition at line 204 of file rsensor.c.

{
       char   linebuf[8192];
       int    nelem = 1000;
       float  *sarr = (float *)malloc(sizeof(float)*nelem);
       FILE   *fp;
       char   *cp;
       int    i;
       
       fp = frlibopen(sfile);
       if (fp == NULL) {
              sprintf(errmsg, "cannot open sensor file '%s'", sfile);
              error(SYSTEM, errmsg);
       }
       fgets(linebuf, sizeof(linebuf), fp);
       if (!strncmp(linebuf, "Elevation ", 10))
              fgets(linebuf, sizeof(linebuf), fp);
                                          /* get phi values */
       sarr[0] = .0f;
       if (strncmp(linebuf, "degrees", 7)) {
              sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
              error(USER, errmsg);
       }
       cp = sskip(linebuf);
       ntp[1] = 0;
       for ( ; ; ) {
              sarr[ntp[1]+1] = atof(cp);
              cp = fskip(cp);
              if (cp == NULL)
                     break;
              ++ntp[1];
       }
       ntp[0] = 0;                        /* get thetas + data */
       while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
              ++ntp[0];
              if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
                     nelem += (nelem>>2) + ntp[1];
                     sarr = (float *)realloc((void *)sarr,
                                   sizeof(float)*nelem);
                     if (sarr == NULL)
                            error(SYSTEM, "out of memory in load_sensor()");
              }
              cp = linebuf;
              i = ntp[0]*(ntp[1]+1);
              for ( ; ; ) {
                     sarr[i] = atof(cp);
                     cp = fskip(cp);
                     if (cp == NULL)
                            break;
                     ++i;
              }
              if (i == ntp[0]*(ntp[1]+1))
                     break;
              if (i != (ntp[0]+1)*(ntp[1]+1)) {
                     sprintf(errmsg,
                     "bad column count near line %d in sensor file '%s'",
                                   ntp[0]+1, sfile);
                     error(USER, errmsg);
              }
       }
       nelem = i;
       fclose(fp);
       errmsg[0] = '\0';                  /* sanity checks */
       if (ntp[0] <= 0)
              sprintf(errmsg, "no data in sensor file '%s'", sfile);
       else if (fabs(sarr[ntp[1]+1]) > FTINY)
              sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
                            sfile);
       else if (fabs(sarr[1]) > FTINY)
              sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
                            sfile);
       else if (sarr[ntp[1]] <= FTINY)
              sprintf(errmsg,
                     "maximum phi must be positive in sensor file '%s'",
                            sfile);
       else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
              sprintf(errmsg,
                     "maximum theta must be positive in sensor file '%s'",
                            sfile);
       if (errmsg[0])
              error(USER, errmsg);
       return((float *)realloc((void *)sarr, sizeof(float)*nelem));
}

Here is the call graph for this function:

Here is the caller graph for this function:

int main ( int  argc,
char *  argv[] 
)

Definition at line 86 of file rsensor.c.

{
       int    doheader = 1;
       int    optwarn = 0;
       int    i, rval;

       progname = argv[0];
                            /* set up rendering defaults */
       rand_samp = 1;
       dstrsrc = 0.65;
       srcsizerat = 0.1;
       directrelay = 3;
       ambounce = 1;
       maxdepth = -10;
                            /* get options from command line */
       for (i = 1; i < argc; i++) {
              while ((rval = expandarg(&argc, &argv, i)) > 0)
                     ;
              if (rval < 0) {
                     sprintf(errmsg, "cannot expand '%s'", argv[i]);
                     error(SYSTEM, errmsg);
              }
              if (argv[i][0] != '-') {
                     if (i >= argc-1)
                            break;        /* final octree argument */
                     if (!ray_pnprocs) {
                            over_options();
                            if (doheader) {      /* print header */
                                   newheader("RADIANCE", stdout);
                                   printargs(argc, argv, stdout);
                                   fputformat("ascii", stdout);
                                   putchar('\n');
                            }
                                          /* start process(es) */
                            ray_pinit(argv[argc-1], nprocs);
                     }
                     comp_sensor(argv[i]);       /* process a sensor file */
                     continue;
              }
              if (argv[i][1] == 'r') {    /* sampling options */
                     if (argv[i][2] == 'd')
                            nsamps = atol(argv[++i]);
                     else if (argv[i][2] == 's')
                            nssamps = atol(argv[++i]);
                     else {
                            sprintf(errmsg, "bad option at '%s'", argv[i]);
                            error(USER, errmsg);
                     }
                     continue;
              }
                                          /* direct component samples */
              if (argv[i][1] == 'd' && argv[i][2] == 'n') {
                     ndsamps = atoi(argv[++i]);
                     continue;
              }
              if (argv[i][1] == 'v') {    /* next sensor view */
                     if (argv[i][2] == 'f') {
                            rval = viewfile(argv[++i], &ourview, NULL);
                            if (rval < 0) {
                                   sprintf(errmsg,
                                   "cannot open view file \"%s\"",
                                                 argv[i]);
                                   error(SYSTEM, errmsg);
                            } else if (rval == 0) {
                                   sprintf(errmsg,
                                          "bad view file \"%s\"",
                                                 argv[i]);
                                   error(USER, errmsg);
                            }
                            continue;
                     }
                     rval = getviewopt(&ourview, argc-i, argv+i);
                     if (rval >= 0) {
                            i += rval;
                            continue;
                     }
                     sprintf(errmsg, "bad view option at '%s'", argv[i]);
                     error(USER, errmsg);
              }
              if (!strcmp(argv[i], "-w")) {      /* toggle warnings */
                     nowarn = !nowarn;
                     continue;
              }
              if (ray_pnprocs) {
                     if (!optwarn++)
                            error(WARNING,
                     "rendering options should appear before first sensor");
              } else if (!strcmp(argv[i], "-defaults")) {
                     print_defaults();
                     return(0);
              }
              if (argv[i][1] == 'h') {    /* header toggle */
                     doheader = !doheader;
                     continue;
              }
              if (!strcmp(argv[i], "-n")) {      /* number of processes */
                     nprocs = atoi(argv[++i]);
                     if (nprocs <= 0)
                            error(USER, "illegal number of processes");
                     continue;
              }
              rval = getrenderopt(argc-i, argv+i);
              if (rval < 0) {
                     sprintf(errmsg, "bad render option at '%s'", argv[i]);
                     error(USER, errmsg);
              }
              i += rval;
       }
       if (!ray_pnprocs)
              error(USER, i<argc ? "missing sensor file" : "missing octree");
       quit(0);
}

Here is the call graph for this function:

static void over_options ( ) [static]

Definition at line 50 of file rsensor.c.

{
       directvis = (ndsamps <= 0);
       do_irrad = 0;
}

Here is the caller graph for this function:

static void print_defaults ( ) [static]

Definition at line 57 of file rsensor.c.

{
       over_options();
       printf("-n %-9d\t\t\t# number of processes\n", nprocs);
       printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
       /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
       printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
       printf("-vp %f %f %f\t# view point\n",
                     ourview.vp[0], ourview.vp[1], ourview.vp[2]);
       printf("-vd %f %f %f\t# view direction\n",
                     ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]);
       printf("-vu %f %f %f\t# view up\n",
                     ourview.vup[0], ourview.vup[1], ourview.vup[2]);
       printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore);
       print_rdefaults();
}

Here is the call graph for this function:

Here is the caller graph for this function:

void quit ( int  ec)

Definition at line 76 of file rsensor.c.

{
       if (ray_pnprocs > 0) /* close children if any */
              ray_pclose(0);              
       exit(ec);
}

Here is the call graph for this function:

static float sens_val ( FVECT  dvec) [static]

Definition at line 434 of file rsensor.c.

{
       FVECT  dv;
       float  theta, phi;
       int    t, p;
       
       dv[2] = DOT(dvec, ourview.vdir);
       theta = (float)((1./DEGREE) * acos(dv[2]));
       if (theta >= maxtheta)
              return(.0f);
       dv[0] = DOT(dvec, ourview.hvec);
       dv[1] = DOT(dvec, ourview.vvec);
       phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
       while (phi < .0f) phi += 360.f;
       t = (int)(theta/maxtheta * sntp[0]);
       p = (int)(phi*(1./360.) * sntp[1]);
                     /* hack for non-uniform sensor grid */
       while (t+1 < sntp[0] && theta >= s_theta(t+1))
              ++t;
       while (t-1 >= 0 && theta <= s_theta(t-1))
              --t;
       while (p+1 < sntp[1] && phi >= s_phi(p+1))
              ++p;
       while (p-1 >= 0 && phi <= s_phi(p-1))
              --p;
       return(s_val(t,p));
}

Here is the caller graph for this function:


Variable Documentation

double gscale = 1.

Definition at line 41 of file rsensor.c.

float maxtheta

Definition at line 36 of file rsensor.c.

int ndsamps = 32

Definition at line 31 of file rsensor.c.

int nowarn

Definition at line 92 of file rcalc.c.

int nphi = 0

Definition at line 40 of file rsensor.c.

int nprocs = 1

Definition at line 32 of file rsensor.c.

unsigned long nsamps = 10000

Definition at line 29 of file rsensor.c.

unsigned long nssamps = 9000

Definition at line 30 of file rsensor.c.

int ntheta = 0

Definition at line 39 of file rsensor.c.

Initial value:
  {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,0.,0.},
                                1.,180.,180.,0.,0.,0.,0.,
                                {0.,0.,0.},{0.,0.,0.},0.,0.}

Definition at line 25 of file rsensor.c.

char* progname

Definition at line 21 of file histo.c.

float* pvals = NULL

Definition at line 38 of file rsensor.c.

const char RCSid[] = "$Id: rsensor.c,v 2.7 2009/12/12 19:01:00 greg Exp $" [static]

Definition at line 2 of file rsensor.c.

float* sensor = NULL

Definition at line 34 of file rsensor.c.

int sntp[2]

Definition at line 35 of file rsensor.c.

float tvals[MAXNT+1]

Definition at line 37 of file rsensor.c.