--- ray/src/util/rsensor.c 2008/02/21 01:22:06 2.1 +++ ray/src/util/rsensor.c 2020/07/24 17:09:33 2.20 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: rsensor.c,v 2.1 2008/02/21 01:22:06 greg Exp $"; +static const char RCSid[] = "$Id: rsensor.c,v 2.20 2020/07/24 17:09:33 greg Exp $"; #endif /* @@ -9,40 +9,55 @@ static const char RCSid[] = "$Id: rsensor.c,v 2.1 2008 */ #include "ray.h" +#include "platform.h" #include "source.h" #include "view.h" #include "random.h" #define DEGREE (PI/180.) -#define MAXNT 180 /* maximum number of theta divisions */ +#define MAXNT 181 /* maximum number of theta divisions */ #define MAXNP 360 /* maximum number of phi divisions */ extern char *progname; /* global argv[0] */ extern int nowarn; /* don't report warnings? */ /* current sensor's perspective */ -VIEW ourview = STDVIEW; +VIEW ourview = {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.}; unsigned long nsamps = 10000; /* desired number of initial samples */ unsigned long nssamps = 9000; /* number of super-samples */ -int ndsamps = 16; /* number of direct samples */ +int ndsamps = 32; /* number of direct samples */ int nprocs = 1; /* number of rendering processes */ float *sensor = NULL; /* current sensor data */ int sntp[2]; /* number of sensor theta and phi angles */ float maxtheta; /* maximum theta value for this sensor */ -float tvals[MAXNT+1]; /* theta values (1-D table of 1-cos(t)) */ -float *pvals = NULL; /* phi values (2-D table in radians) */ +float tvals[MAXNT+1]; /* theta prob. values (1-D table of 1-cos(t)) */ +float *pvals = NULL; /* phi prob. values (2-D table in radians) */ int ntheta = 0; /* polar angle divisions */ int nphi = 0; /* azimuthal angle divisions */ double gscale = 1.; /* global scaling value */ +#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)] + static void comp_sensor(char *sfile); static void -print_defaults() +over_options() /* overriding options */ { + directvis = (ndsamps <= 0); + do_irrad = 0; +} + +static void +print_defaults() /* print out default parameters */ +{ + 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); */ @@ -57,6 +72,19 @@ print_defaults() print_rdefaults(); } + +void +quit(ec) /* make sure exit is called */ +int ec; +{ + if (ray_pnprocs > 0) /* close children if any */ + ray_pclose(0); + else if (ray_pnprocs < 0) + _exit(ec); /* avoid flush in child */ + exit(ec); +} + + int main( int argc, @@ -64,43 +92,41 @@ main( ) { int doheader = 1; + int optwarn = 0; int i, rval; progname = argv[0]; /* set up rendering defaults */ - dstrsrc = 0.25; + rand_samp = 1; + dstrsrc = 0.65; + srcsizerat = 0.1; directrelay = 3; ambounce = 1; - /* just asking defaults? */ - if (argc == 2 && !strcmp(argv[1], "-defaults")) { - print_defaults(); - return(0); - } - /* check octree */ - if (argc < 2 || argv[argc-1][0] == '-') - error(USER, "missing octree argument"); + maxdepth = -10; /* get options from command line */ - for (i = 1; i < argc-1; i++) { + 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] != '-') { /* process a sensor file */ + if (argv[i][0] != '-') { + if (i >= argc-1) + break; /* final octree argument */ if (!ray_pnprocs) { - /* overriding options */ - directvis = (ndsamps <= 0); - do_irrad = 0; + 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); + if (strcmp(argv[argc-1], ".")) + ray_pinit(argv[argc-1], nprocs); } - comp_sensor(argv[i]); + comp_sensor(argv[i]); /* process a sensor file */ continue; } if (argv[i][1] == 'r') { /* sampling options */ @@ -143,12 +169,13 @@ main( sprintf(errmsg, "bad view option at '%s'", argv[i]); error(USER, errmsg); } - if (!strcmp(argv[i], "-w")) { /* turn off warnings */ - nowarn = 1; + if (!strcmp(argv[i], "-w")) { /* toggle warnings */ + nowarn = !nowarn; continue; } if (ray_pnprocs) { - error(WARNING, + if (!optwarn++) + error(WARNING, "rendering options should appear before first sensor"); } else if (!strcmp(argv[i], "-defaults")) { print_defaults(); @@ -171,6 +198,8 @@ main( } i += rval; } + if (sensor == NULL) + error(USER, i 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]+FTINY) { + sprintf(errmsg, + "Phi values not monotinically increasing in sensor file '%s'", + sfile); + error(USER, errmsg); + } ++ntp[1]; } + warnedneg = 0; ntp[0] = 0; /* get thetas + data */ while (fgets(linebuf, sizeof(linebuf), fp) != NULL) { ++ntp[0]; @@ -228,10 +266,26 @@ load_sensor( cp = fskip(cp); if (cp == NULL) break; + if (sarr[i] < .0) { + if (!warnedneg++) { + sprintf(errmsg, + "Negative value(s) in sensor file '%s' (ignored)\n", sfile); + error(WARNING, errmsg); + } + sarr[i] = .0; + } else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1)) + last_pos_val = i; ++i; } - if (i == ntp[0]*(ntp[1]+1)) + if (i == ntp[0]*(ntp[1]+1)) /* empty line? */ break; + if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <= + sarr[(ntp[0]-1)*(ntp[1]+1)]) { + sprintf(errmsg, + "Theta values not monotinically increasing in sensor file '%s'", + sfile); + error(USER, errmsg); + } if (i != (ntp[0]+1)*(ntp[1]+1)) { sprintf(errmsg, "bad column count near line %d in sensor file '%s'", @@ -239,24 +293,26 @@ load_sensor( error(USER, errmsg); } } - nelem = i; + /* truncate zero region */ + ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1; + nelem = (ntp[0]+1)*(ntp[1]+1); fclose(fp); errmsg[0] = '\0'; /* sanity checks */ - if (ntp[0] <= 0) - sprintf(errmsg, "no data in sensor file '%s'", sfile); + if (!last_pos_val) + sprintf(errmsg, "no positive sensor values in 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) + else if (sarr[ntp[1]] < 270.-FTINY) sprintf(errmsg, - "maximum phi must be positive in sensor file '%s'", + "maximum phi must be 270 or greater in sensor file '%s'", sfile); - else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY) + else if (sarr[ntp[1]] >= 360.-FTINY) sprintf(errmsg, - "maximum theta must be positive in sensor file '%s'", + "maximum phi must be less than 360 in sensor file '%s'", sfile); if (errmsg[0]) error(USER, errmsg); @@ -282,6 +338,8 @@ init_ptable( if (pvals != NULL) free((void *)pvals); if (sfile == NULL || !*sfile) { + sensor = NULL; + sntp[0] = sntp[1] = 0; pvals = NULL; ntheta = nphi = 0; return; @@ -299,41 +357,44 @@ init_ptable( error(INTERNAL, errmsg); } /* compute boundary angles */ - maxtheta = 1.5f*sensor[sntp[0]*(sntp[1]+1)] - - 0.5f*sensor[sntp[0]*sntp[1]]; + maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2)); + if (maxtheta > PI) + maxtheta = PI; thdiv[0] = .0; for (t = 1; t < sntp[0]; t++) - thdiv[t] = DEGREE/2.*(sensor[t*(sntp[1]+1)] + - sensor[(t+1)*(sntp[1]+1)]); - thdiv[sntp[0]] = maxtheta*DEGREE; - phdiv[0] = .0; + thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t)); + thdiv[sntp[0]] = maxtheta; + phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1)); for (p = 1; p < sntp[1]; p++) - phdiv[p] = DEGREE/2.*(sensor[p] + sensor[p+1]); - phdiv[sntp[1]] = 2.*PI; + phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p)); + phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2)); /* size our table */ - tsize = 1. - cos(maxtheta*DEGREE); + tsize = 1. - cos(maxtheta); psize = PI*tsize/maxtheta; if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */ samptot = sntp[0]*sntp[1]; - ntheta = (int)(sqrt(samptot*tsize/psize) + 0.5); + ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1; if (ntheta > MAXNT) ntheta = MAXNT; nphi = samptot/ntheta; - pvals = (float *)malloc(sizeof(float)*ntheta*(nphi+1)); + pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(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 = sensor + (i+1)*(sntp[1]+1) + 1; - rowsum[i] = 0.; + rowp = &s_val(i,0); + rowsum[i] = 1e-20; 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]; } - tvals[0] = .0f; - for (i = 1; i < ntheta; i++) { + if (gscale <= FTINY) { + sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile); + error(USER, errmsg); + } + 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) @@ -343,26 +404,37 @@ init_ptable( frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale); tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) + frac*cos(thdiv[t+1]) ); - pvals[i*(nphi+1)] = .0f; + /* offset b/c sensor values are centered */ + if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) { + frac -= 0.5; + } else { + frac += 0.5; + --t; + } + pvals[i*(nphi+1)] = phdiv[0]; for (j = 1; j < nphi; j++) { prob = (double)j / (double)nphi; - rowp = sensor + t*(sntp[1]+1) + 1; - rowp1 = rowp + sntp[1]+1; - for (p = 0; p < sntp[1]; p++) { - if ((prob -= (1.-frac)*rowp[p]/rowsum[t-1] + - frac*rowp1[p]/rowsum[t]) <= .0) + 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-1] - + frac*rowp1[p]/rowsum[t]); - pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] + - frac1*phdiv[p+1]; + if (p >= sntp[1]) { /* should never happen? */ + p = sntp[1] - 1; + prob = .5; } + 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); + pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]]; } + /* duplicate final row */ + memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1), + sizeof(*pvals)*(nphi+1)); + tvals[0] = .0f; tvals[ntheta] = (float)tsize; } @@ -408,27 +480,41 @@ sens_val( int t, p; dv[2] = DOT(dvec, ourview.vdir); - theta = (float)((1./DEGREE) * acos(dv[2])); + theta = 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; + phi = atan2(-dv[0], dv[1]); + while (phi < .0f) phi += (float)(2.*PI); t = (int)(theta/maxtheta * sntp[0]); - p = (int)(phi*(1./360.) * sntp[1]); + p = (int)(phi*(1./(2.*PI)) * sntp[1]); /* hack for non-uniform sensor grid */ - while (t+1 < sntp[0] && theta >= sensor[(t+2)*(sntp[1]+1)]) + theta *= (float)(1./DEGREE); + phi *= (float)(1./DEGREE); + while (t+1 < sntp[0] && theta >= s_theta(t+1)) ++t; - while (t-1 >= 0 && theta < sensor[t*(sntp[1]+1)]) + while (t-1 >= 0 && theta <= s_theta(t-1)) --t; - while (p+1 < sntp[1] && phi >= sensor[p+2]) + while (p+1 < sntp[1] && phi >= s_phi(p+1)) ++p; - while (p-1 >= 0 && phi < sensor[p]) + while (p-1 >= 0 && phi <= s_phi(p-1)) --p; - return(sensor[t*(sntp[1]+1) + p + 1]); + return(s_val(t,p)); } +/* Print origin and direction */ +static void +print_ray( + FVECT rorg, + FVECT rdir +) +{ + printf("%.6g %.6g %.6g %.8f %.8f %.8f\n", + rorg[0], rorg[1], rorg[2], + rdir[0], rdir[1], rdir[2]); +} + /* Compute sensor output */ static void comp_sensor( @@ -441,6 +527,7 @@ comp_sensor( int nt, np; COLOR vsum; RAY rr; + double sf; int i, j; /* set view */ ourview.type = VT_ANG; @@ -451,41 +538,64 @@ comp_sensor( error(USER, err); /* assign probability table */ init_ptable(sfile); - /* do Monte Carlo sampling */ + /* stratified MC sampling */ setcolor(vsum, .0f, .0f, .0f); nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5); np = nsamps/nt; - VCOPY(rr.rorg, ourview.vp); - rr.rmax = .0; + sf = gscale/nsamps; for (i = 0; i < nt; i++) - for (j =0; j < np; j++) { - get_direc(rr.rdir, (i+frandom())/nt, - (j + frandom())/np); - rayorigin(&rr, PRIMARY, NULL, NULL); + 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); + if (!ray_pnprocs) { + print_ray(rr.rorg, rr.rdir); + continue; + } + rr.rmax = .0; + rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL); + scalecolor(rr.rcoef, sf); if (ray_pqueue(&rr) == 1) addcolor(vsum, rr.rcol); } - /* finish MC calculation */ - while (ray_presult(&rr, 0) > 0) - addcolor(vsum, rr.rcol); - scalecolor(vsum, gscale/(nt*np)); - /* compute direct component */ + /* 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); + if (!ray_pnprocs) { + print_ray(rr.rorg, rr.rdir); + continue; + } + rr.rmax = .0; + rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL); + scalecolor(rr.rcoef, sf); + if (ray_pqueue(&rr) == 1) + addcolor(vsum, rr.rcol); + } + if (!ray_pnprocs) /* just printing rays */ + return; + /* scale partial result */ + scalecolor(vsum, sf); + /* add direct component */ for (i = ndirs; i-- > 0; ) { SRCINDEX si; initsrcindex(&si); while (srcray(&rr, NULL, &si)) { - double d = sens_val(rr.rdir); - if (d <= FTINY) + sf = sens_val(rr.rdir); + if (sf <= FTINY) continue; - d *= si.dom/ndirs; - scalecolor(rr.rcoef, d); + sf *= si.dom/ndirs; + scalecolor(rr.rcoef, sf); if (ray_pqueue(&rr) == 1) { multcolor(rr.rcol, rr.rcoef); addcolor(vsum, rr.rcol); } } } - /* finish direct calculation */ + /* finish our calculation */ while (ray_presult(&rr, 0) > 0) { multcolor(rr.rcol, rr.rcoef); addcolor(vsum, rr.rcol);