--- ray/src/util/rsensor.c 2010/09/26 15:41:46 2.8 +++ ray/src/util/rsensor.c 2022/03/15 18:05:03 2.21 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: rsensor.c,v 2.8 2010/09/26 15:41:46 greg Exp $"; +static const char RCSid[] = "$Id: rsensor.c,v 2.21 2022/03/15 18:05:03 greg Exp $"; #endif /* @@ -9,13 +9,14 @@ static const char RCSid[] = "$Id: rsensor.c,v 2.8 2010 */ #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] */ @@ -26,16 +27,15 @@ VIEW ourview = {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,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 */ +long nsamps = 10000; /* desired number of initial 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 */ @@ -59,7 +59,6 @@ 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); */ 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]); @@ -77,7 +76,9 @@ quit(ec) /* make sure exit is called */ int ec; { if (ray_pnprocs > 0) /* close children if any */ - ray_pclose(0); + ray_pclose(0); + else if (ray_pnprocs < 0) + _exit(ec); /* avoid flush in child */ exit(ec); } @@ -129,8 +130,6 @@ main( 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); @@ -207,7 +206,9 @@ load_sensor( char *sfile ) { + int warnedneg; char linebuf[8192]; + int last_pos_val = 0; int nelem = 1000; float *sarr = (float *)malloc(sizeof(float)*nelem); FILE *fp; @@ -235,8 +236,15 @@ load_sensor( cp = fskip(cp); if (cp == NULL) break; + if (ntp[1] > 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]; @@ -254,10 +262,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'", @@ -265,24 +289,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); @@ -295,7 +321,7 @@ init_ptable( char *sfile ) { - int samptot = nsamps; + long samptot = nsamps; float *rowp, *rowp1; double rowsum[MAXNT], rowomega[MAXNT]; double thdiv[MAXNT+1], phdiv[MAXNP+1]; @@ -327,37 +353,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); + 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 = &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++) @@ -369,32 +401,35 @@ 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) + if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) { frac -= 0.5; - else if (t >= sntp[0]-1 || frac < 0.5) { + } 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]); - 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; } @@ -441,16 +476,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)) @@ -513,7 +550,7 @@ comp_sensor( continue; } rr.rmax = .0; - rayorigin(&rr, PRIMARY, NULL, NULL); + rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL); scalecolor(rr.rcoef, sf); if (ray_pqueue(&rr) == 1) addcolor(vsum, rr.rcol); @@ -529,7 +566,7 @@ comp_sensor( continue; } rr.rmax = .0; - rayorigin(&rr, PRIMARY, NULL, NULL); + rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL); scalecolor(rr.rcoef, sf); if (ray_pqueue(&rr) == 1) addcolor(vsum, rr.rcol);