--- ray/src/util/rsensor.c 2008/02/22 18:04:02 2.3 +++ ray/src/util/rsensor.c 2013/10/16 04:35:50 2.13 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: rsensor.c,v 2.3 2008/02/22 18:04:02 greg Exp $"; +static const char RCSid[] = "$Id: rsensor.c,v 2.13 2013/10/16 04:35:50 greg Exp $"; #endif /* @@ -22,7 +22,9 @@ 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 */ @@ -32,8 +34,8 @@ 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 */ @@ -69,6 +71,17 @@ print_defaults() /* print out default parameters */ print_rdefaults(); } + +void +quit(ec) /* make sure exit is called */ +int ec; +{ + if (ray_pnprocs > 0) /* close children if any */ + ray_pclose(0); + exit(ec); +} + + int main( int argc, @@ -82,7 +95,7 @@ main( progname = argv[0]; /* set up rendering defaults */ rand_samp = 1; - dstrsrc = 0.5; + dstrsrc = 0.65; srcsizerat = 0.1; directrelay = 3; ambounce = 1; @@ -101,12 +114,14 @@ main( 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); + if (strcmp(argv[argc-1], ".")) + ray_pinit(argv[argc-1], nprocs); } comp_sensor(argv[i]); /* process a sensor file */ continue; @@ -180,7 +195,7 @@ main( } i += rval; } - if (!ray_pnprocs) + if (sensor == NULL) error(USER, i 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]) { + 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]; @@ -239,10 +262,25 @@ load_sensor( cp = fskip(cp); if (cp == NULL) break; + if (i && sarr[i] < .0) { + if (!warnedneg++) { + sprintf(errmsg, + "Negative value(s) in sensor file '%s' (ignored)\n", sfile); + error(WARNING, errmsg); + } + sarr[i] = .0; + } ++i; } if (i == ntp[0]*(ntp[1]+1)) 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'", @@ -312,37 +350,43 @@ init_ptable( error(INTERNAL, errmsg); } /* compute boundary angles */ - maxtheta = 1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2); + 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.*(s_theta(t-1) + s_theta(t)); - thdiv[sntp[0]] = maxtheta*DEGREE; - phdiv[0] = .0; + 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.*(s_phi(p-1) + s_phi(p)); - phdiv[sntp[1]] = 2.*PI; + 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); - psize = PI*tsize/(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((double)samptot*tsize/psize) + 0.5); 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 = &s_val(i,0); - rowsum[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]; } + 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++) @@ -354,37 +398,32 @@ init_ptable( 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) { + if ((t < sntp[0]-1) & (frac >= 0.5)) { + if ((frac -= 0.5) < 0) + frac = 0; + } else { frac += 0.5; --t; } - pvals[i*(nphi+1)] = .0f; + pvals[i*(nphi+1)] = phdiv[0]; 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++) { + 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]); - if (p <= 0 || frac1 > 0.5) - frac1 -= 0.5; - else if (p >= sntp[1]-1 || frac1 < 0.5) { - frac1 += 0.5; - --p; - } - 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]]; } tvals[0] = .0f; tvals[ntheta] = (float)tsize; @@ -432,16 +471,18 @@ 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 */ + theta *= (float)(1./DEGREE); + phi *= (float)(1./DEGREE); while (t+1 < sntp[0] && theta >= s_theta(t+1)) ++t; while (t-1 >= 0 && theta <= s_theta(t-1)) @@ -453,6 +494,18 @@ sens_val( 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( @@ -465,6 +518,7 @@ comp_sensor( int nt, np; COLOR vsum; RAY rr; + double sf; int i, j; /* set view */ ourview.type = VT_ANG; @@ -475,40 +529,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++) { + VCOPY(rr.rorg, ourview.vp); get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np); - rayorigin(&rr, PRIMARY, NULL, NULL); + 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);