Back to index

radiance  4R0+20100331
rsensor.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: rsensor.c,v 2.7 2009/12/12 19:01:00 greg Exp $";
00003 #endif
00004 
00005 /*
00006  * Compute sensor signal based on spatial sensitivity.
00007  *
00008  *     Created Feb 2008 for Architectural Energy Corp.
00009  */
00010 
00011 #include "ray.h"
00012 #include "source.h"
00013 #include "view.h"
00014 #include "random.h"
00015 
00016 #define DEGREE              (PI/180.)
00017 
00018 #define       MAXNT         180    /* maximum number of theta divisions */
00019 #define MAXNP        360    /* maximum number of phi divisions */
00020 
00021 extern char   *progname;    /* global argv[0] */
00022 extern int    nowarn;         /* don't report warnings? */
00023 
00024                             /* current sensor's perspective */
00025 VIEW          ourview =  {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,0.,0.},
00026                                 1.,180.,180.,0.,0.,0.,0.,
00027                                 {0.,0.,0.},{0.,0.,0.},0.,0.};
00028 
00029 unsigned long nsamps = 10000;      /* desired number of initial samples */
00030 unsigned long nssamps = 9000;      /* number of super-samples */
00031 int           ndsamps = 32; /* number of direct samples */
00032 int           nprocs = 1;     /* number of rendering processes */
00033 
00034 float         *sensor = NULL;      /* current sensor data */
00035 int           sntp[2];      /* number of sensor theta and phi angles */
00036 float         maxtheta;     /* maximum theta value for this sensor */
00037 float         tvals[MAXNT+1];      /* theta values (1-D table of 1-cos(t)) */
00038 float         *pvals = NULL;       /* phi values (2-D table in radians) */
00039 int           ntheta = 0;   /* polar angle divisions */
00040 int           nphi = 0;     /* azimuthal angle divisions */
00041 double        gscale = 1.;  /* global scaling value */
00042 
00043 #define       s_theta(t)    sensor[(t+1)*(sntp[1]+1)]
00044 #define s_phi(p)     sensor[(p)+1]
00045 #define s_val(t,p)   sensor[(p)+1+(t+1)*(sntp[1]+1)]
00046 
00047 static void   comp_sensor(char *sfile);
00048 
00049 static void
00050 over_options()                     /* overriding options */
00051 {
00052        directvis = (ndsamps <= 0);
00053        do_irrad = 0;
00054 }
00055 
00056 static void
00057 print_defaults()            /* print out default parameters */
00058 {
00059        over_options();
00060        printf("-n %-9d\t\t\t# number of processes\n", nprocs);
00061        printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
00062        /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
00063        printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
00064        printf("-vp %f %f %f\t# view point\n",
00065                      ourview.vp[0], ourview.vp[1], ourview.vp[2]);
00066        printf("-vd %f %f %f\t# view direction\n",
00067                      ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]);
00068        printf("-vu %f %f %f\t# view up\n",
00069                      ourview.vup[0], ourview.vup[1], ourview.vup[2]);
00070        printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore);
00071        print_rdefaults();
00072 }
00073 
00074 
00075 void
00076 quit(ec)                    /* make sure exit is called */
00077 int    ec;
00078 {
00079        if (ray_pnprocs > 0) /* close children if any */
00080               ray_pclose(0);              
00081        exit(ec);
00082 }
00083 
00084 
00085 int
00086 main(
00087        int    argc,
00088        char   *argv[]
00089 )
00090 {
00091        int    doheader = 1;
00092        int    optwarn = 0;
00093        int    i, rval;
00094 
00095        progname = argv[0];
00096                             /* set up rendering defaults */
00097        rand_samp = 1;
00098        dstrsrc = 0.65;
00099        srcsizerat = 0.1;
00100        directrelay = 3;
00101        ambounce = 1;
00102        maxdepth = -10;
00103                             /* get options from command line */
00104        for (i = 1; i < argc; i++) {
00105               while ((rval = expandarg(&argc, &argv, i)) > 0)
00106                      ;
00107               if (rval < 0) {
00108                      sprintf(errmsg, "cannot expand '%s'", argv[i]);
00109                      error(SYSTEM, errmsg);
00110               }
00111               if (argv[i][0] != '-') {
00112                      if (i >= argc-1)
00113                             break;        /* final octree argument */
00114                      if (!ray_pnprocs) {
00115                             over_options();
00116                             if (doheader) {      /* print header */
00117                                    newheader("RADIANCE", stdout);
00118                                    printargs(argc, argv, stdout);
00119                                    fputformat("ascii", stdout);
00120                                    putchar('\n');
00121                             }
00122                                           /* start process(es) */
00123                             ray_pinit(argv[argc-1], nprocs);
00124                      }
00125                      comp_sensor(argv[i]);       /* process a sensor file */
00126                      continue;
00127               }
00128               if (argv[i][1] == 'r') {    /* sampling options */
00129                      if (argv[i][2] == 'd')
00130                             nsamps = atol(argv[++i]);
00131                      else if (argv[i][2] == 's')
00132                             nssamps = atol(argv[++i]);
00133                      else {
00134                             sprintf(errmsg, "bad option at '%s'", argv[i]);
00135                             error(USER, errmsg);
00136                      }
00137                      continue;
00138               }
00139                                           /* direct component samples */
00140               if (argv[i][1] == 'd' && argv[i][2] == 'n') {
00141                      ndsamps = atoi(argv[++i]);
00142                      continue;
00143               }
00144               if (argv[i][1] == 'v') {    /* next sensor view */
00145                      if (argv[i][2] == 'f') {
00146                             rval = viewfile(argv[++i], &ourview, NULL);
00147                             if (rval < 0) {
00148                                    sprintf(errmsg,
00149                                    "cannot open view file \"%s\"",
00150                                                  argv[i]);
00151                                    error(SYSTEM, errmsg);
00152                             } else if (rval == 0) {
00153                                    sprintf(errmsg,
00154                                           "bad view file \"%s\"",
00155                                                  argv[i]);
00156                                    error(USER, errmsg);
00157                             }
00158                             continue;
00159                      }
00160                      rval = getviewopt(&ourview, argc-i, argv+i);
00161                      if (rval >= 0) {
00162                             i += rval;
00163                             continue;
00164                      }
00165                      sprintf(errmsg, "bad view option at '%s'", argv[i]);
00166                      error(USER, errmsg);
00167               }
00168               if (!strcmp(argv[i], "-w")) {      /* toggle warnings */
00169                      nowarn = !nowarn;
00170                      continue;
00171               }
00172               if (ray_pnprocs) {
00173                      if (!optwarn++)
00174                             error(WARNING,
00175                      "rendering options should appear before first sensor");
00176               } else if (!strcmp(argv[i], "-defaults")) {
00177                      print_defaults();
00178                      return(0);
00179               }
00180               if (argv[i][1] == 'h') {    /* header toggle */
00181                      doheader = !doheader;
00182                      continue;
00183               }
00184               if (!strcmp(argv[i], "-n")) {      /* number of processes */
00185                      nprocs = atoi(argv[++i]);
00186                      if (nprocs <= 0)
00187                             error(USER, "illegal number of processes");
00188                      continue;
00189               }
00190               rval = getrenderopt(argc-i, argv+i);
00191               if (rval < 0) {
00192                      sprintf(errmsg, "bad render option at '%s'", argv[i]);
00193                      error(USER, errmsg);
00194               }
00195               i += rval;
00196        }
00197        if (!ray_pnprocs)
00198               error(USER, i<argc ? "missing sensor file" : "missing octree");
00199        quit(0);
00200 }
00201 
00202 /* Load sensor sensitivities (first row and column are angles) */
00203 static float *
00204 load_sensor(
00205        int    ntp[2],
00206        char   *sfile
00207 )
00208 {
00209        char   linebuf[8192];
00210        int    nelem = 1000;
00211        float  *sarr = (float *)malloc(sizeof(float)*nelem);
00212        FILE   *fp;
00213        char   *cp;
00214        int    i;
00215        
00216        fp = frlibopen(sfile);
00217        if (fp == NULL) {
00218               sprintf(errmsg, "cannot open sensor file '%s'", sfile);
00219               error(SYSTEM, errmsg);
00220        }
00221        fgets(linebuf, sizeof(linebuf), fp);
00222        if (!strncmp(linebuf, "Elevation ", 10))
00223               fgets(linebuf, sizeof(linebuf), fp);
00224                                           /* get phi values */
00225        sarr[0] = .0f;
00226        if (strncmp(linebuf, "degrees", 7)) {
00227               sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
00228               error(USER, errmsg);
00229        }
00230        cp = sskip(linebuf);
00231        ntp[1] = 0;
00232        for ( ; ; ) {
00233               sarr[ntp[1]+1] = atof(cp);
00234               cp = fskip(cp);
00235               if (cp == NULL)
00236                      break;
00237               ++ntp[1];
00238        }
00239        ntp[0] = 0;                        /* get thetas + data */
00240        while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
00241               ++ntp[0];
00242               if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
00243                      nelem += (nelem>>2) + ntp[1];
00244                      sarr = (float *)realloc((void *)sarr,
00245                                    sizeof(float)*nelem);
00246                      if (sarr == NULL)
00247                             error(SYSTEM, "out of memory in load_sensor()");
00248               }
00249               cp = linebuf;
00250               i = ntp[0]*(ntp[1]+1);
00251               for ( ; ; ) {
00252                      sarr[i] = atof(cp);
00253                      cp = fskip(cp);
00254                      if (cp == NULL)
00255                             break;
00256                      ++i;
00257               }
00258               if (i == ntp[0]*(ntp[1]+1))
00259                      break;
00260               if (i != (ntp[0]+1)*(ntp[1]+1)) {
00261                      sprintf(errmsg,
00262                      "bad column count near line %d in sensor file '%s'",
00263                                    ntp[0]+1, sfile);
00264                      error(USER, errmsg);
00265               }
00266        }
00267        nelem = i;
00268        fclose(fp);
00269        errmsg[0] = '\0';                  /* sanity checks */
00270        if (ntp[0] <= 0)
00271               sprintf(errmsg, "no data in sensor file '%s'", sfile);
00272        else if (fabs(sarr[ntp[1]+1]) > FTINY)
00273               sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
00274                             sfile);
00275        else if (fabs(sarr[1]) > FTINY)
00276               sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
00277                             sfile);
00278        else if (sarr[ntp[1]] <= FTINY)
00279               sprintf(errmsg,
00280                      "maximum phi must be positive in sensor file '%s'",
00281                             sfile);
00282        else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
00283               sprintf(errmsg,
00284                      "maximum theta must be positive in sensor file '%s'",
00285                             sfile);
00286        if (errmsg[0])
00287               error(USER, errmsg);
00288        return((float *)realloc((void *)sarr, sizeof(float)*nelem));
00289 }
00290 
00291 /* Initialize probability table */
00292 static void
00293 init_ptable(
00294        char   *sfile
00295 )
00296 {
00297        int    samptot = nsamps;
00298        float  *rowp, *rowp1;
00299        double rowsum[MAXNT], rowomega[MAXNT];
00300        double thdiv[MAXNT+1], phdiv[MAXNP+1];
00301        double tsize, psize;
00302        double prob, frac, frac1;
00303        int    i, j, t, p;
00304                                    /* free old table */
00305        if (sensor != NULL)
00306               free((void *)sensor);
00307        if (pvals != NULL)
00308               free((void *)pvals);
00309        if (sfile == NULL || !*sfile) {
00310               sensor = NULL;
00311               sntp[0] = sntp[1] = 0;
00312               pvals = NULL;
00313               ntheta = nphi = 0;
00314               return;
00315        }
00316                                    /* load sensor table */
00317        sensor = load_sensor(sntp, sfile);
00318        if (sntp[0] > MAXNT) {
00319               sprintf(errmsg, "Too many theta rows in sensor file '%s'",
00320                             sfile);
00321               error(INTERNAL, errmsg);
00322        }
00323        if (sntp[1] > MAXNP) {
00324               sprintf(errmsg, "Too many phi columns in sensor file '%s'",
00325                             sfile);
00326               error(INTERNAL, errmsg);
00327        }
00328                                    /* compute boundary angles */
00329        maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2);
00330        thdiv[0] = .0;
00331        for (t = 1; t < sntp[0]; t++)
00332               thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
00333        thdiv[sntp[0]] = maxtheta*DEGREE;
00334        phdiv[0] = .0;
00335        for (p = 1; p < sntp[1]; p++)
00336               phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
00337        phdiv[sntp[1]] = 2.*PI;
00338                                    /* size our table */
00339        tsize = 1. - cos(maxtheta*DEGREE);
00340        psize = PI*tsize/(maxtheta*DEGREE);
00341        if (sntp[0]*sntp[1] < samptot)     /* don't overdo resolution */
00342               samptot = sntp[0]*sntp[1];
00343        ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
00344        if (ntheta > MAXNT)
00345               ntheta = MAXNT;
00346        nphi = samptot/ntheta;
00347        pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1));
00348        if (pvals == NULL)
00349               error(SYSTEM, "out of memory in init_ptable()");
00350        gscale = .0;                /* compute our inverse table */
00351        for (i = 0; i < sntp[0]; i++) {
00352               rowp = &s_val(i,0);
00353               rowsum[i] = 0.;
00354               for (j = 0; j < sntp[1]; j++)
00355                      rowsum[i] += *rowp++;
00356               rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
00357               rowomega[i] *= 2.*PI / (double)sntp[1];
00358               gscale += rowsum[i] * rowomega[i];
00359        }
00360        for (i = 0; i < ntheta; i++) {
00361               prob = (double)i / (double)ntheta;
00362               for (t = 0; t < sntp[0]; t++)
00363                      if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
00364                             break;
00365               if (t >= sntp[0])
00366                      error(INTERNAL, "code error 1 in init_ptable()");
00367               frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
00368               tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
00369                                           frac*cos(thdiv[t+1]) );
00370                             /* offset b/c sensor values are centered */
00371               if (t <= 0 || frac > 0.5)
00372                      frac -= 0.5;
00373               else if (t >= sntp[0]-1 || frac < 0.5) {
00374                      frac += 0.5;
00375                      --t;
00376               }
00377               pvals[i*(nphi+1)] = .0f;
00378               for (j = 1; j < nphi; j++) {
00379                      prob = (double)j / (double)nphi;
00380                      rowp = &s_val(t,0);
00381                      rowp1 = &s_val(t+1,0);
00382                      for (p = 0; p < sntp[1]; p++) {
00383                             if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
00384                                        frac*rowp1[p]/rowsum[t+1]) <= .0)
00385                                    break;
00386                             if (p >= sntp[1])
00387                                    error(INTERNAL,
00388                                        "code error 2 in init_ptable()");
00389                             frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
00390                                           + frac*rowp1[p]/rowsum[t+1]);
00391                             pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
00392                                                  frac1*phdiv[p+1];
00393                      }
00394               }
00395               pvals[i*(nphi+1) + nphi] = (float)(2.*PI);
00396        }
00397        tvals[0] = .0f;
00398        tvals[ntheta] = (float)tsize;
00399 }
00400 
00401 /* Get normalized direction from random variables in [0,1) range */
00402 static void
00403 get_direc(
00404        FVECT dvec,
00405        double x,
00406        double y
00407 )
00408 {
00409        double xfrac = x*ntheta;
00410        int    tndx = (int)xfrac;
00411        double yfrac = y*nphi;
00412        int    pndx = (int)yfrac;
00413        double rad, phi;
00414        FVECT  dv;
00415        int    i;
00416 
00417        xfrac -= (double)tndx;
00418        yfrac -= (double)pndx;
00419        pndx += tndx*(nphi+1);
00420        
00421        dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]);
00422        rad = sqrt(1. - dv[2]*dv[2]);
00423        phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1];
00424        dv[0] = -rad*sin(phi);
00425        dv[1] = rad*cos(phi);
00426        for (i = 3; i--; )
00427               dvec[i] = dv[0]*ourview.hvec[i] +
00428                             dv[1]*ourview.vvec[i] +
00429                             dv[2]*ourview.vdir[i] ;
00430 }
00431 
00432 /* Get sensor value in the specified direction (normalized) */
00433 static float
00434 sens_val(
00435        FVECT  dvec
00436 )
00437 {
00438        FVECT  dv;
00439        float  theta, phi;
00440        int    t, p;
00441        
00442        dv[2] = DOT(dvec, ourview.vdir);
00443        theta = (float)((1./DEGREE) * acos(dv[2]));
00444        if (theta >= maxtheta)
00445               return(.0f);
00446        dv[0] = DOT(dvec, ourview.hvec);
00447        dv[1] = DOT(dvec, ourview.vvec);
00448        phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
00449        while (phi < .0f) phi += 360.f;
00450        t = (int)(theta/maxtheta * sntp[0]);
00451        p = (int)(phi*(1./360.) * sntp[1]);
00452                      /* hack for non-uniform sensor grid */
00453        while (t+1 < sntp[0] && theta >= s_theta(t+1))
00454               ++t;
00455        while (t-1 >= 0 && theta <= s_theta(t-1))
00456               --t;
00457        while (p+1 < sntp[1] && phi >= s_phi(p+1))
00458               ++p;
00459        while (p-1 >= 0 && phi <= s_phi(p-1))
00460               --p;
00461        return(s_val(t,p));
00462 }
00463 
00464 /* Compute sensor output */
00465 static void
00466 comp_sensor(
00467        char *sfile
00468 )
00469 {
00470        int    ndirs = dstrsrc > FTINY ? ndsamps :
00471                             ndsamps > 0 ? 1 : 0;
00472        char   *err;
00473        int    nt, np;
00474        COLOR  vsum;
00475        RAY    rr;
00476        double sf;
00477        int    i, j;
00478                                           /* set view */
00479        ourview.type = VT_ANG;
00480        ourview.horiz = ourview.vert = 180.;
00481        ourview.hoff = ourview.voff = .0;
00482        err = setview(&ourview);
00483        if (err != NULL)
00484               error(USER, err);
00485                                           /* assign probability table */
00486        init_ptable(sfile);
00487                                           /* stratified MC sampling */
00488        setcolor(vsum, .0f, .0f, .0f);
00489        nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
00490        np = nsamps/nt;
00491        sf = gscale/nsamps;
00492        for (i = 0; i < nt; i++)
00493               for (j = 0; j < np; j++) {
00494                      VCOPY(rr.rorg, ourview.vp);
00495                      get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
00496                      if (ourview.vfore > FTINY)
00497                             VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
00498                      rr.rmax = .0;
00499                      rayorigin(&rr, PRIMARY, NULL, NULL);
00500                      scalecolor(rr.rcoef, sf);
00501                      if (ray_pqueue(&rr) == 1)
00502                             addcolor(vsum, rr.rcol);
00503               }
00504                                           /* remaining rays pure MC */
00505        for (i = nsamps - nt*np; i-- > 0; ) {
00506               VCOPY(rr.rorg, ourview.vp);
00507               get_direc(rr.rdir, frandom(), frandom());
00508               if (ourview.vfore > FTINY)
00509                      VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
00510               rr.rmax = .0;
00511               rayorigin(&rr, PRIMARY, NULL, NULL);
00512               scalecolor(rr.rcoef, sf);
00513               if (ray_pqueue(&rr) == 1)
00514                      addcolor(vsum, rr.rcol);
00515        }
00516                                           /* scale partial result */
00517        scalecolor(vsum, sf);
00518                                           /* add direct component */
00519        for (i = ndirs; i-- > 0; ) {
00520               SRCINDEX      si;
00521               initsrcindex(&si);
00522               while (srcray(&rr, NULL, &si)) {
00523                      sf = sens_val(rr.rdir);
00524                      if (sf <= FTINY)
00525                             continue;
00526                      sf *= si.dom/ndirs;
00527                      scalecolor(rr.rcoef, sf);
00528                      if (ray_pqueue(&rr) == 1) {
00529                             multcolor(rr.rcol, rr.rcoef);
00530                             addcolor(vsum, rr.rcol);
00531                      }
00532               }
00533        }
00534                                           /* finish our calculation */
00535        while (ray_presult(&rr, 0) > 0) {
00536               multcolor(rr.rcol, rr.rcoef);
00537               addcolor(vsum, rr.rcol);
00538        }
00539                                           /* print our result */
00540        printf("%.4e %.4e %.4e\n", colval(vsum,RED),
00541                             colval(vsum,GRN), colval(vsum,BLU));
00542 }