| 1 |
greg |
2.1 |
#ifndef lint |
| 2 |
greg |
2.24 |
static const char RCSid[] = "$Id: rsensor.c,v 2.23 2023/12/16 03:30:35 greg Exp $"; |
| 3 |
greg |
2.1 |
#endif |
| 4 |
|
|
|
| 5 |
|
|
/* |
| 6 |
|
|
* Compute sensor signal based on spatial sensitivity. |
| 7 |
|
|
* |
| 8 |
|
|
* Created Feb 2008 for Architectural Energy Corp. |
| 9 |
|
|
*/ |
| 10 |
|
|
|
| 11 |
|
|
#include "ray.h" |
| 12 |
greg |
2.20 |
#include "platform.h" |
| 13 |
greg |
2.1 |
#include "source.h" |
| 14 |
greg |
2.24 |
#include "func.h" |
| 15 |
greg |
2.1 |
#include "view.h" |
| 16 |
|
|
#include "random.h" |
| 17 |
|
|
|
| 18 |
|
|
#define DEGREE (PI/180.) |
| 19 |
|
|
|
| 20 |
greg |
2.15 |
#define MAXNT 181 /* maximum number of theta divisions */ |
| 21 |
greg |
2.1 |
#define MAXNP 360 /* maximum number of phi divisions */ |
| 22 |
|
|
|
| 23 |
|
|
extern char *progname; /* global argv[0] */ |
| 24 |
|
|
extern int nowarn; /* don't report warnings? */ |
| 25 |
|
|
|
| 26 |
|
|
/* current sensor's perspective */ |
| 27 |
greg |
2.5 |
VIEW ourview = {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,0.,0.}, |
| 28 |
|
|
1.,180.,180.,0.,0.,0.,0., |
| 29 |
|
|
{0.,0.,0.},{0.,0.,0.},0.,0.}; |
| 30 |
greg |
2.1 |
|
| 31 |
greg |
2.21 |
long nsamps = 10000; /* desired number of initial samples */ |
| 32 |
greg |
2.2 |
int ndsamps = 32; /* number of direct samples */ |
| 33 |
greg |
2.1 |
int nprocs = 1; /* number of rendering processes */ |
| 34 |
|
|
|
| 35 |
|
|
float *sensor = NULL; /* current sensor data */ |
| 36 |
|
|
int sntp[2]; /* number of sensor theta and phi angles */ |
| 37 |
|
|
float maxtheta; /* maximum theta value for this sensor */ |
| 38 |
greg |
2.9 |
float tvals[MAXNT+1]; /* theta prob. values (1-D table of 1-cos(t)) */ |
| 39 |
|
|
float *pvals = NULL; /* phi prob. values (2-D table in radians) */ |
| 40 |
greg |
2.1 |
int ntheta = 0; /* polar angle divisions */ |
| 41 |
|
|
int nphi = 0; /* azimuthal angle divisions */ |
| 42 |
|
|
double gscale = 1.; /* global scaling value */ |
| 43 |
|
|
|
| 44 |
greg |
2.2 |
#define s_theta(t) sensor[(t+1)*(sntp[1]+1)] |
| 45 |
|
|
#define s_phi(p) sensor[(p)+1] |
| 46 |
|
|
#define s_val(t,p) sensor[(p)+1+(t+1)*(sntp[1]+1)] |
| 47 |
|
|
|
| 48 |
greg |
2.1 |
static void comp_sensor(char *sfile); |
| 49 |
|
|
|
| 50 |
|
|
static void |
| 51 |
greg |
2.3 |
over_options() /* overriding options */ |
| 52 |
greg |
2.1 |
{ |
| 53 |
greg |
2.3 |
directvis = (ndsamps <= 0); |
| 54 |
|
|
do_irrad = 0; |
| 55 |
|
|
} |
| 56 |
|
|
|
| 57 |
|
|
static void |
| 58 |
|
|
print_defaults() /* print out default parameters */ |
| 59 |
|
|
{ |
| 60 |
|
|
over_options(); |
| 61 |
greg |
2.1 |
printf("-n %-9d\t\t\t# number of processes\n", nprocs); |
| 62 |
|
|
printf("-rd %-9ld\t\t\t# ray directions\n", nsamps); |
| 63 |
|
|
printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps); |
| 64 |
|
|
printf("-vp %f %f %f\t# view point\n", |
| 65 |
|
|
ourview.vp[0], ourview.vp[1], ourview.vp[2]); |
| 66 |
|
|
printf("-vd %f %f %f\t# view direction\n", |
| 67 |
|
|
ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]); |
| 68 |
|
|
printf("-vu %f %f %f\t# view up\n", |
| 69 |
|
|
ourview.vup[0], ourview.vup[1], ourview.vup[2]); |
| 70 |
|
|
printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore); |
| 71 |
|
|
print_rdefaults(); |
| 72 |
|
|
} |
| 73 |
|
|
|
| 74 |
greg |
2.7 |
|
| 75 |
|
|
void |
| 76 |
|
|
quit(ec) /* make sure exit is called */ |
| 77 |
|
|
int ec; |
| 78 |
|
|
{ |
| 79 |
|
|
if (ray_pnprocs > 0) /* close children if any */ |
| 80 |
greg |
2.19 |
ray_pclose(0); |
| 81 |
|
|
else if (ray_pnprocs < 0) |
| 82 |
|
|
_exit(ec); /* avoid flush in child */ |
| 83 |
greg |
2.7 |
exit(ec); |
| 84 |
|
|
} |
| 85 |
|
|
|
| 86 |
|
|
|
| 87 |
greg |
2.1 |
int |
| 88 |
|
|
main( |
| 89 |
|
|
int argc, |
| 90 |
|
|
char *argv[] |
| 91 |
|
|
) |
| 92 |
|
|
{ |
| 93 |
|
|
int doheader = 1; |
| 94 |
greg |
2.3 |
int optwarn = 0; |
| 95 |
greg |
2.1 |
int i, rval; |
| 96 |
|
|
|
| 97 |
|
|
progname = argv[0]; |
| 98 |
|
|
/* set up rendering defaults */ |
| 99 |
greg |
2.2 |
rand_samp = 1; |
| 100 |
greg |
2.6 |
dstrsrc = 0.65; |
| 101 |
greg |
2.2 |
srcsizerat = 0.1; |
| 102 |
greg |
2.1 |
directrelay = 3; |
| 103 |
|
|
ambounce = 1; |
| 104 |
greg |
2.2 |
maxdepth = -10; |
| 105 |
greg |
2.24 |
/* initialize calcomp routines */ |
| 106 |
|
|
initfunc(); |
| 107 |
greg |
2.1 |
/* get options from command line */ |
| 108 |
greg |
2.3 |
for (i = 1; i < argc; i++) { |
| 109 |
greg |
2.1 |
while ((rval = expandarg(&argc, &argv, i)) > 0) |
| 110 |
|
|
; |
| 111 |
|
|
if (rval < 0) { |
| 112 |
|
|
sprintf(errmsg, "cannot expand '%s'", argv[i]); |
| 113 |
|
|
error(SYSTEM, errmsg); |
| 114 |
|
|
} |
| 115 |
greg |
2.3 |
if (argv[i][0] != '-') { |
| 116 |
|
|
if (i >= argc-1) |
| 117 |
|
|
break; /* final octree argument */ |
| 118 |
greg |
2.1 |
if (!ray_pnprocs) { |
| 119 |
greg |
2.3 |
over_options(); |
| 120 |
greg |
2.1 |
if (doheader) { /* print header */ |
| 121 |
greg |
2.6 |
newheader("RADIANCE", stdout); |
| 122 |
greg |
2.1 |
printargs(argc, argv, stdout); |
| 123 |
|
|
fputformat("ascii", stdout); |
| 124 |
|
|
putchar('\n'); |
| 125 |
|
|
} |
| 126 |
|
|
/* start process(es) */ |
| 127 |
greg |
2.8 |
if (strcmp(argv[argc-1], ".")) |
| 128 |
|
|
ray_pinit(argv[argc-1], nprocs); |
| 129 |
greg |
2.1 |
} |
| 130 |
greg |
2.3 |
comp_sensor(argv[i]); /* process a sensor file */ |
| 131 |
greg |
2.1 |
continue; |
| 132 |
|
|
} |
| 133 |
|
|
if (argv[i][1] == 'r') { /* sampling options */ |
| 134 |
|
|
if (argv[i][2] == 'd') |
| 135 |
|
|
nsamps = atol(argv[++i]); |
| 136 |
|
|
else { |
| 137 |
|
|
sprintf(errmsg, "bad option at '%s'", argv[i]); |
| 138 |
|
|
error(USER, errmsg); |
| 139 |
|
|
} |
| 140 |
|
|
continue; |
| 141 |
|
|
} |
| 142 |
|
|
/* direct component samples */ |
| 143 |
|
|
if (argv[i][1] == 'd' && argv[i][2] == 'n') { |
| 144 |
|
|
ndsamps = atoi(argv[++i]); |
| 145 |
|
|
continue; |
| 146 |
|
|
} |
| 147 |
|
|
if (argv[i][1] == 'v') { /* next sensor view */ |
| 148 |
|
|
if (argv[i][2] == 'f') { |
| 149 |
|
|
rval = viewfile(argv[++i], &ourview, NULL); |
| 150 |
|
|
if (rval < 0) { |
| 151 |
|
|
sprintf(errmsg, |
| 152 |
|
|
"cannot open view file \"%s\"", |
| 153 |
|
|
argv[i]); |
| 154 |
|
|
error(SYSTEM, errmsg); |
| 155 |
|
|
} else if (rval == 0) { |
| 156 |
|
|
sprintf(errmsg, |
| 157 |
|
|
"bad view file \"%s\"", |
| 158 |
|
|
argv[i]); |
| 159 |
|
|
error(USER, errmsg); |
| 160 |
|
|
} |
| 161 |
|
|
continue; |
| 162 |
|
|
} |
| 163 |
|
|
rval = getviewopt(&ourview, argc-i, argv+i); |
| 164 |
|
|
if (rval >= 0) { |
| 165 |
|
|
i += rval; |
| 166 |
|
|
continue; |
| 167 |
|
|
} |
| 168 |
|
|
sprintf(errmsg, "bad view option at '%s'", argv[i]); |
| 169 |
|
|
error(USER, errmsg); |
| 170 |
|
|
} |
| 171 |
greg |
2.3 |
if (!strcmp(argv[i], "-w")) { /* toggle warnings */ |
| 172 |
|
|
nowarn = !nowarn; |
| 173 |
greg |
2.1 |
continue; |
| 174 |
|
|
} |
| 175 |
|
|
if (ray_pnprocs) { |
| 176 |
greg |
2.3 |
if (!optwarn++) |
| 177 |
|
|
error(WARNING, |
| 178 |
greg |
2.1 |
"rendering options should appear before first sensor"); |
| 179 |
|
|
} else if (!strcmp(argv[i], "-defaults")) { |
| 180 |
|
|
print_defaults(); |
| 181 |
|
|
return(0); |
| 182 |
|
|
} |
| 183 |
|
|
if (argv[i][1] == 'h') { /* header toggle */ |
| 184 |
|
|
doheader = !doheader; |
| 185 |
|
|
continue; |
| 186 |
|
|
} |
| 187 |
|
|
if (!strcmp(argv[i], "-n")) { /* number of processes */ |
| 188 |
|
|
nprocs = atoi(argv[++i]); |
| 189 |
|
|
if (nprocs <= 0) |
| 190 |
|
|
error(USER, "illegal number of processes"); |
| 191 |
|
|
continue; |
| 192 |
|
|
} |
| 193 |
|
|
rval = getrenderopt(argc-i, argv+i); |
| 194 |
|
|
if (rval < 0) { |
| 195 |
|
|
sprintf(errmsg, "bad render option at '%s'", argv[i]); |
| 196 |
|
|
error(USER, errmsg); |
| 197 |
|
|
} |
| 198 |
|
|
i += rval; |
| 199 |
|
|
} |
| 200 |
greg |
2.8 |
if (sensor == NULL) |
| 201 |
greg |
2.3 |
error(USER, i<argc ? "missing sensor file" : "missing octree"); |
| 202 |
greg |
2.1 |
quit(0); |
| 203 |
|
|
} |
| 204 |
|
|
|
| 205 |
|
|
/* Load sensor sensitivities (first row and column are angles) */ |
| 206 |
|
|
static float * |
| 207 |
|
|
load_sensor( |
| 208 |
|
|
int ntp[2], |
| 209 |
|
|
char *sfile |
| 210 |
|
|
) |
| 211 |
|
|
{ |
| 212 |
greg |
2.10 |
int warnedneg; |
| 213 |
greg |
2.1 |
char linebuf[8192]; |
| 214 |
greg |
2.16 |
int last_pos_val = 0; |
| 215 |
greg |
2.1 |
int nelem = 1000; |
| 216 |
|
|
float *sarr = (float *)malloc(sizeof(float)*nelem); |
| 217 |
|
|
FILE *fp; |
| 218 |
|
|
char *cp; |
| 219 |
|
|
int i; |
| 220 |
|
|
|
| 221 |
|
|
fp = frlibopen(sfile); |
| 222 |
|
|
if (fp == NULL) { |
| 223 |
|
|
sprintf(errmsg, "cannot open sensor file '%s'", sfile); |
| 224 |
|
|
error(SYSTEM, errmsg); |
| 225 |
|
|
} |
| 226 |
|
|
fgets(linebuf, sizeof(linebuf), fp); |
| 227 |
|
|
if (!strncmp(linebuf, "Elevation ", 10)) |
| 228 |
|
|
fgets(linebuf, sizeof(linebuf), fp); |
| 229 |
|
|
/* get phi values */ |
| 230 |
|
|
sarr[0] = .0f; |
| 231 |
|
|
if (strncmp(linebuf, "degrees", 7)) { |
| 232 |
|
|
sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile); |
| 233 |
|
|
error(USER, errmsg); |
| 234 |
|
|
} |
| 235 |
|
|
cp = sskip(linebuf); |
| 236 |
|
|
ntp[1] = 0; |
| 237 |
|
|
for ( ; ; ) { |
| 238 |
|
|
sarr[ntp[1]+1] = atof(cp); |
| 239 |
|
|
cp = fskip(cp); |
| 240 |
|
|
if (cp == NULL) |
| 241 |
|
|
break; |
| 242 |
greg |
2.18 |
if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]+FTINY) { |
| 243 |
greg |
2.9 |
sprintf(errmsg, |
| 244 |
|
|
"Phi values not monotinically increasing in sensor file '%s'", |
| 245 |
|
|
sfile); |
| 246 |
|
|
error(USER, errmsg); |
| 247 |
|
|
} |
| 248 |
greg |
2.1 |
++ntp[1]; |
| 249 |
|
|
} |
| 250 |
greg |
2.10 |
warnedneg = 0; |
| 251 |
greg |
2.1 |
ntp[0] = 0; /* get thetas + data */ |
| 252 |
|
|
while (fgets(linebuf, sizeof(linebuf), fp) != NULL) { |
| 253 |
|
|
++ntp[0]; |
| 254 |
|
|
if ((ntp[0]+1)*(ntp[1]+1) > nelem) { |
| 255 |
|
|
nelem += (nelem>>2) + ntp[1]; |
| 256 |
|
|
sarr = (float *)realloc((void *)sarr, |
| 257 |
|
|
sizeof(float)*nelem); |
| 258 |
|
|
if (sarr == NULL) |
| 259 |
|
|
error(SYSTEM, "out of memory in load_sensor()"); |
| 260 |
|
|
} |
| 261 |
|
|
cp = linebuf; |
| 262 |
|
|
i = ntp[0]*(ntp[1]+1); |
| 263 |
|
|
for ( ; ; ) { |
| 264 |
|
|
sarr[i] = atof(cp); |
| 265 |
|
|
cp = fskip(cp); |
| 266 |
|
|
if (cp == NULL) |
| 267 |
|
|
break; |
| 268 |
greg |
2.17 |
if (sarr[i] < .0) { |
| 269 |
greg |
2.10 |
if (!warnedneg++) { |
| 270 |
|
|
sprintf(errmsg, |
| 271 |
|
|
"Negative value(s) in sensor file '%s' (ignored)\n", sfile); |
| 272 |
|
|
error(WARNING, errmsg); |
| 273 |
|
|
} |
| 274 |
|
|
sarr[i] = .0; |
| 275 |
greg |
2.17 |
} else if (sarr[i] > FTINY && i > ntp[0]*(ntp[1]+1)) |
| 276 |
greg |
2.16 |
last_pos_val = i; |
| 277 |
greg |
2.1 |
++i; |
| 278 |
|
|
} |
| 279 |
greg |
2.17 |
if (i == ntp[0]*(ntp[1]+1)) /* empty line? */ |
| 280 |
greg |
2.1 |
break; |
| 281 |
greg |
2.9 |
if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <= |
| 282 |
|
|
sarr[(ntp[0]-1)*(ntp[1]+1)]) { |
| 283 |
|
|
sprintf(errmsg, |
| 284 |
|
|
"Theta values not monotinically increasing in sensor file '%s'", |
| 285 |
|
|
sfile); |
| 286 |
|
|
error(USER, errmsg); |
| 287 |
|
|
} |
| 288 |
greg |
2.1 |
if (i != (ntp[0]+1)*(ntp[1]+1)) { |
| 289 |
|
|
sprintf(errmsg, |
| 290 |
|
|
"bad column count near line %d in sensor file '%s'", |
| 291 |
|
|
ntp[0]+1, sfile); |
| 292 |
|
|
error(USER, errmsg); |
| 293 |
|
|
} |
| 294 |
|
|
} |
| 295 |
greg |
2.16 |
/* truncate zero region */ |
| 296 |
|
|
ntp[0] = (last_pos_val + ntp[1])/(ntp[1]+1) - 1; |
| 297 |
|
|
nelem = (ntp[0]+1)*(ntp[1]+1); |
| 298 |
greg |
2.1 |
fclose(fp); |
| 299 |
|
|
errmsg[0] = '\0'; /* sanity checks */ |
| 300 |
greg |
2.16 |
if (!last_pos_val) |
| 301 |
|
|
sprintf(errmsg, "no positive sensor values in file '%s'", sfile); |
| 302 |
greg |
2.1 |
else if (fabs(sarr[ntp[1]+1]) > FTINY) |
| 303 |
|
|
sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'", |
| 304 |
|
|
sfile); |
| 305 |
|
|
else if (fabs(sarr[1]) > FTINY) |
| 306 |
|
|
sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'", |
| 307 |
|
|
sfile); |
| 308 |
greg |
2.18 |
else if (sarr[ntp[1]] < 270.-FTINY) |
| 309 |
greg |
2.1 |
sprintf(errmsg, |
| 310 |
greg |
2.18 |
"maximum phi must be 270 or greater in sensor file '%s'", |
| 311 |
|
|
sfile); |
| 312 |
|
|
else if (sarr[ntp[1]] >= 360.-FTINY) |
| 313 |
|
|
sprintf(errmsg, |
| 314 |
|
|
"maximum phi must be less than 360 in sensor file '%s'", |
| 315 |
greg |
2.1 |
sfile); |
| 316 |
|
|
if (errmsg[0]) |
| 317 |
|
|
error(USER, errmsg); |
| 318 |
|
|
return((float *)realloc((void *)sarr, sizeof(float)*nelem)); |
| 319 |
|
|
} |
| 320 |
|
|
|
| 321 |
|
|
/* Initialize probability table */ |
| 322 |
|
|
static void |
| 323 |
|
|
init_ptable( |
| 324 |
|
|
char *sfile |
| 325 |
|
|
) |
| 326 |
|
|
{ |
| 327 |
greg |
2.21 |
long samptot = nsamps; |
| 328 |
greg |
2.1 |
float *rowp, *rowp1; |
| 329 |
|
|
double rowsum[MAXNT], rowomega[MAXNT]; |
| 330 |
|
|
double thdiv[MAXNT+1], phdiv[MAXNP+1]; |
| 331 |
|
|
double tsize, psize; |
| 332 |
|
|
double prob, frac, frac1; |
| 333 |
|
|
int i, j, t, p; |
| 334 |
|
|
/* free old table */ |
| 335 |
|
|
if (sensor != NULL) |
| 336 |
|
|
free((void *)sensor); |
| 337 |
|
|
if (pvals != NULL) |
| 338 |
|
|
free((void *)pvals); |
| 339 |
|
|
if (sfile == NULL || !*sfile) { |
| 340 |
greg |
2.2 |
sensor = NULL; |
| 341 |
|
|
sntp[0] = sntp[1] = 0; |
| 342 |
greg |
2.1 |
pvals = NULL; |
| 343 |
|
|
ntheta = nphi = 0; |
| 344 |
|
|
return; |
| 345 |
|
|
} |
| 346 |
|
|
/* load sensor table */ |
| 347 |
|
|
sensor = load_sensor(sntp, sfile); |
| 348 |
|
|
if (sntp[0] > MAXNT) { |
| 349 |
|
|
sprintf(errmsg, "Too many theta rows in sensor file '%s'", |
| 350 |
|
|
sfile); |
| 351 |
|
|
error(INTERNAL, errmsg); |
| 352 |
|
|
} |
| 353 |
|
|
if (sntp[1] > MAXNP) { |
| 354 |
|
|
sprintf(errmsg, "Too many phi columns in sensor file '%s'", |
| 355 |
|
|
sfile); |
| 356 |
|
|
error(INTERNAL, errmsg); |
| 357 |
|
|
} |
| 358 |
|
|
/* compute boundary angles */ |
| 359 |
greg |
2.9 |
maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2)); |
| 360 |
|
|
if (maxtheta > PI) |
| 361 |
|
|
maxtheta = PI; |
| 362 |
greg |
2.1 |
thdiv[0] = .0; |
| 363 |
|
|
for (t = 1; t < sntp[0]; t++) |
| 364 |
greg |
2.2 |
thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t)); |
| 365 |
greg |
2.9 |
thdiv[sntp[0]] = maxtheta; |
| 366 |
|
|
phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1)); |
| 367 |
greg |
2.1 |
for (p = 1; p < sntp[1]; p++) |
| 368 |
greg |
2.2 |
phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p)); |
| 369 |
greg |
2.9 |
phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2)); |
| 370 |
greg |
2.1 |
/* size our table */ |
| 371 |
greg |
2.9 |
tsize = 1. - cos(maxtheta); |
| 372 |
|
|
psize = PI*tsize/maxtheta; |
| 373 |
greg |
2.1 |
if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */ |
| 374 |
|
|
samptot = sntp[0]*sntp[1]; |
| 375 |
greg |
2.17 |
ntheta = (int)(sqrt((double)samptot*tsize/psize)*sntp[0]/sntp[1]) + 1; |
| 376 |
greg |
2.1 |
if (ntheta > MAXNT) |
| 377 |
|
|
ntheta = MAXNT; |
| 378 |
|
|
nphi = samptot/ntheta; |
| 379 |
greg |
2.9 |
pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1)); |
| 380 |
greg |
2.1 |
if (pvals == NULL) |
| 381 |
|
|
error(SYSTEM, "out of memory in init_ptable()"); |
| 382 |
|
|
gscale = .0; /* compute our inverse table */ |
| 383 |
|
|
for (i = 0; i < sntp[0]; i++) { |
| 384 |
greg |
2.2 |
rowp = &s_val(i,0); |
| 385 |
greg |
2.9 |
rowsum[i] = 1e-20; |
| 386 |
greg |
2.1 |
for (j = 0; j < sntp[1]; j++) |
| 387 |
|
|
rowsum[i] += *rowp++; |
| 388 |
|
|
rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]); |
| 389 |
|
|
rowomega[i] *= 2.*PI / (double)sntp[1]; |
| 390 |
|
|
gscale += rowsum[i] * rowomega[i]; |
| 391 |
|
|
} |
| 392 |
greg |
2.9 |
if (gscale <= FTINY) { |
| 393 |
|
|
sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile); |
| 394 |
|
|
error(USER, errmsg); |
| 395 |
|
|
} |
| 396 |
greg |
2.2 |
for (i = 0; i < ntheta; i++) { |
| 397 |
greg |
2.1 |
prob = (double)i / (double)ntheta; |
| 398 |
|
|
for (t = 0; t < sntp[0]; t++) |
| 399 |
|
|
if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0) |
| 400 |
|
|
break; |
| 401 |
|
|
if (t >= sntp[0]) |
| 402 |
|
|
error(INTERNAL, "code error 1 in init_ptable()"); |
| 403 |
|
|
frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale); |
| 404 |
|
|
tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) + |
| 405 |
|
|
frac*cos(thdiv[t+1]) ); |
| 406 |
greg |
2.2 |
/* offset b/c sensor values are centered */ |
| 407 |
greg |
2.14 |
if ((t < sntp[0]-1) & (!t | (frac >= 0.5))) { |
| 408 |
|
|
frac -= 0.5; |
| 409 |
greg |
2.13 |
} else { |
| 410 |
greg |
2.2 |
frac += 0.5; |
| 411 |
|
|
--t; |
| 412 |
|
|
} |
| 413 |
greg |
2.9 |
pvals[i*(nphi+1)] = phdiv[0]; |
| 414 |
greg |
2.1 |
for (j = 1; j < nphi; j++) { |
| 415 |
|
|
prob = (double)j / (double)nphi; |
| 416 |
greg |
2.2 |
rowp = &s_val(t,0); |
| 417 |
|
|
rowp1 = &s_val(t+1,0); |
| 418 |
greg |
2.9 |
for (p = 0; p < sntp[1]; p++) |
| 419 |
greg |
2.2 |
if ((prob -= (1.-frac)*rowp[p]/rowsum[t] + |
| 420 |
|
|
frac*rowp1[p]/rowsum[t+1]) <= .0) |
| 421 |
greg |
2.1 |
break; |
| 422 |
greg |
2.10 |
if (p >= sntp[1]) { /* should never happen? */ |
| 423 |
greg |
2.9 |
p = sntp[1] - 1; |
| 424 |
|
|
prob = .5; |
| 425 |
greg |
2.1 |
} |
| 426 |
greg |
2.9 |
frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t] |
| 427 |
|
|
+ frac*rowp1[p]/rowsum[t+1]); |
| 428 |
|
|
pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] + |
| 429 |
|
|
frac1*phdiv[p+1]; |
| 430 |
greg |
2.1 |
} |
| 431 |
greg |
2.9 |
pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]]; |
| 432 |
greg |
2.1 |
} |
| 433 |
greg |
2.15 |
/* duplicate final row */ |
| 434 |
|
|
memcpy(pvals+ntheta*(nphi+1), pvals+(ntheta-1)*(nphi+1), |
| 435 |
|
|
sizeof(*pvals)*(nphi+1)); |
| 436 |
greg |
2.2 |
tvals[0] = .0f; |
| 437 |
greg |
2.1 |
tvals[ntheta] = (float)tsize; |
| 438 |
|
|
} |
| 439 |
|
|
|
| 440 |
|
|
/* Get normalized direction from random variables in [0,1) range */ |
| 441 |
|
|
static void |
| 442 |
|
|
get_direc( |
| 443 |
|
|
FVECT dvec, |
| 444 |
|
|
double x, |
| 445 |
|
|
double y |
| 446 |
|
|
) |
| 447 |
|
|
{ |
| 448 |
|
|
double xfrac = x*ntheta; |
| 449 |
|
|
int tndx = (int)xfrac; |
| 450 |
|
|
double yfrac = y*nphi; |
| 451 |
|
|
int pndx = (int)yfrac; |
| 452 |
|
|
double rad, phi; |
| 453 |
|
|
FVECT dv; |
| 454 |
|
|
int i; |
| 455 |
|
|
|
| 456 |
|
|
xfrac -= (double)tndx; |
| 457 |
|
|
yfrac -= (double)pndx; |
| 458 |
|
|
pndx += tndx*(nphi+1); |
| 459 |
|
|
|
| 460 |
|
|
dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]); |
| 461 |
|
|
rad = sqrt(1. - dv[2]*dv[2]); |
| 462 |
|
|
phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1]; |
| 463 |
|
|
dv[0] = -rad*sin(phi); |
| 464 |
|
|
dv[1] = rad*cos(phi); |
| 465 |
|
|
for (i = 3; i--; ) |
| 466 |
|
|
dvec[i] = dv[0]*ourview.hvec[i] + |
| 467 |
|
|
dv[1]*ourview.vvec[i] + |
| 468 |
|
|
dv[2]*ourview.vdir[i] ; |
| 469 |
|
|
} |
| 470 |
|
|
|
| 471 |
|
|
/* Get sensor value in the specified direction (normalized) */ |
| 472 |
|
|
static float |
| 473 |
|
|
sens_val( |
| 474 |
|
|
FVECT dvec |
| 475 |
|
|
) |
| 476 |
|
|
{ |
| 477 |
|
|
FVECT dv; |
| 478 |
|
|
float theta, phi; |
| 479 |
|
|
int t, p; |
| 480 |
|
|
|
| 481 |
|
|
dv[2] = DOT(dvec, ourview.vdir); |
| 482 |
greg |
2.11 |
theta = acos(dv[2]); |
| 483 |
greg |
2.1 |
if (theta >= maxtheta) |
| 484 |
|
|
return(.0f); |
| 485 |
|
|
dv[0] = DOT(dvec, ourview.hvec); |
| 486 |
|
|
dv[1] = DOT(dvec, ourview.vvec); |
| 487 |
greg |
2.11 |
phi = atan2(-dv[0], dv[1]); |
| 488 |
|
|
while (phi < .0f) phi += (float)(2.*PI); |
| 489 |
greg |
2.1 |
t = (int)(theta/maxtheta * sntp[0]); |
| 490 |
greg |
2.11 |
p = (int)(phi*(1./(2.*PI)) * sntp[1]); |
| 491 |
greg |
2.1 |
/* hack for non-uniform sensor grid */ |
| 492 |
greg |
2.11 |
theta *= (float)(1./DEGREE); |
| 493 |
|
|
phi *= (float)(1./DEGREE); |
| 494 |
greg |
2.2 |
while (t+1 < sntp[0] && theta >= s_theta(t+1)) |
| 495 |
greg |
2.1 |
++t; |
| 496 |
greg |
2.2 |
while (t-1 >= 0 && theta <= s_theta(t-1)) |
| 497 |
greg |
2.1 |
--t; |
| 498 |
greg |
2.2 |
while (p+1 < sntp[1] && phi >= s_phi(p+1)) |
| 499 |
greg |
2.1 |
++p; |
| 500 |
greg |
2.2 |
while (p-1 >= 0 && phi <= s_phi(p-1)) |
| 501 |
greg |
2.1 |
--p; |
| 502 |
greg |
2.2 |
return(s_val(t,p)); |
| 503 |
greg |
2.1 |
} |
| 504 |
|
|
|
| 505 |
greg |
2.8 |
/* Print origin and direction */ |
| 506 |
|
|
static void |
| 507 |
|
|
print_ray( |
| 508 |
|
|
FVECT rorg, |
| 509 |
|
|
FVECT rdir |
| 510 |
|
|
) |
| 511 |
|
|
{ |
| 512 |
|
|
printf("%.6g %.6g %.6g %.8f %.8f %.8f\n", |
| 513 |
|
|
rorg[0], rorg[1], rorg[2], |
| 514 |
|
|
rdir[0], rdir[1], rdir[2]); |
| 515 |
|
|
} |
| 516 |
|
|
|
| 517 |
greg |
2.1 |
/* Compute sensor output */ |
| 518 |
|
|
static void |
| 519 |
|
|
comp_sensor( |
| 520 |
|
|
char *sfile |
| 521 |
|
|
) |
| 522 |
|
|
{ |
| 523 |
|
|
int ndirs = dstrsrc > FTINY ? ndsamps : |
| 524 |
|
|
ndsamps > 0 ? 1 : 0; |
| 525 |
|
|
char *err; |
| 526 |
|
|
int nt, np; |
| 527 |
greg |
2.22 |
SCOLOR vsum; |
| 528 |
greg |
2.1 |
RAY rr; |
| 529 |
greg |
2.6 |
double sf; |
| 530 |
greg |
2.1 |
int i, j; |
| 531 |
|
|
/* set view */ |
| 532 |
|
|
ourview.type = VT_ANG; |
| 533 |
|
|
ourview.horiz = ourview.vert = 180.; |
| 534 |
|
|
ourview.hoff = ourview.voff = .0; |
| 535 |
|
|
err = setview(&ourview); |
| 536 |
|
|
if (err != NULL) |
| 537 |
|
|
error(USER, err); |
| 538 |
|
|
/* assign probability table */ |
| 539 |
|
|
init_ptable(sfile); |
| 540 |
greg |
2.6 |
/* stratified MC sampling */ |
| 541 |
greg |
2.22 |
scolorblack(vsum); |
| 542 |
greg |
2.1 |
nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5); |
| 543 |
|
|
np = nsamps/nt; |
| 544 |
greg |
2.6 |
sf = gscale/nsamps; |
| 545 |
greg |
2.1 |
for (i = 0; i < nt; i++) |
| 546 |
greg |
2.2 |
for (j = 0; j < np; j++) { |
| 547 |
greg |
2.5 |
VCOPY(rr.rorg, ourview.vp); |
| 548 |
greg |
2.2 |
get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np); |
| 549 |
greg |
2.5 |
if (ourview.vfore > FTINY) |
| 550 |
|
|
VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore); |
| 551 |
greg |
2.8 |
if (!ray_pnprocs) { |
| 552 |
|
|
print_ray(rr.rorg, rr.rdir); |
| 553 |
|
|
continue; |
| 554 |
|
|
} |
| 555 |
greg |
2.5 |
rr.rmax = .0; |
| 556 |
greg |
2.12 |
rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL); |
| 557 |
greg |
2.22 |
scalescolor(rr.rcoef, sf); |
| 558 |
greg |
2.1 |
if (ray_pqueue(&rr) == 1) |
| 559 |
greg |
2.22 |
saddscolor(vsum, rr.rcol); |
| 560 |
greg |
2.1 |
} |
| 561 |
greg |
2.6 |
/* remaining rays pure MC */ |
| 562 |
|
|
for (i = nsamps - nt*np; i-- > 0; ) { |
| 563 |
|
|
VCOPY(rr.rorg, ourview.vp); |
| 564 |
|
|
get_direc(rr.rdir, frandom(), frandom()); |
| 565 |
|
|
if (ourview.vfore > FTINY) |
| 566 |
|
|
VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore); |
| 567 |
greg |
2.8 |
if (!ray_pnprocs) { |
| 568 |
|
|
print_ray(rr.rorg, rr.rdir); |
| 569 |
|
|
continue; |
| 570 |
|
|
} |
| 571 |
greg |
2.6 |
rr.rmax = .0; |
| 572 |
greg |
2.12 |
rayorigin(&rr, PRIMARY|SPECULAR, NULL, NULL); |
| 573 |
greg |
2.22 |
scalescolor(rr.rcoef, sf); |
| 574 |
greg |
2.6 |
if (ray_pqueue(&rr) == 1) |
| 575 |
greg |
2.22 |
saddscolor(vsum, rr.rcol); |
| 576 |
greg |
2.6 |
} |
| 577 |
greg |
2.8 |
if (!ray_pnprocs) /* just printing rays */ |
| 578 |
|
|
return; |
| 579 |
greg |
2.6 |
/* scale partial result */ |
| 580 |
greg |
2.22 |
scalescolor(vsum, sf); |
| 581 |
greg |
2.6 |
/* add direct component */ |
| 582 |
greg |
2.1 |
for (i = ndirs; i-- > 0; ) { |
| 583 |
|
|
SRCINDEX si; |
| 584 |
|
|
initsrcindex(&si); |
| 585 |
|
|
while (srcray(&rr, NULL, &si)) { |
| 586 |
greg |
2.6 |
sf = sens_val(rr.rdir); |
| 587 |
|
|
if (sf <= FTINY) |
| 588 |
greg |
2.1 |
continue; |
| 589 |
greg |
2.6 |
sf *= si.dom/ndirs; |
| 590 |
greg |
2.22 |
scalescolor(rr.rcoef, sf); |
| 591 |
greg |
2.1 |
if (ray_pqueue(&rr) == 1) { |
| 592 |
greg |
2.22 |
smultscolor(rr.rcol, rr.rcoef); |
| 593 |
|
|
saddscolor(vsum, rr.rcol); |
| 594 |
greg |
2.1 |
} |
| 595 |
|
|
} |
| 596 |
|
|
} |
| 597 |
greg |
2.23 |
while (ray_presult(&rr, 0) > 0) { /* finish our calculation */ |
| 598 |
greg |
2.22 |
smultscolor(rr.rcol, rr.rcoef); |
| 599 |
|
|
saddscolor(vsum, rr.rcol); |
| 600 |
greg |
2.1 |
} |
| 601 |
greg |
2.23 |
for (i = 0; i < NCSAMP; i++) /* print our result */ |
| 602 |
|
|
printf(" %.4e", vsum[i]); |
| 603 |
|
|
fputc('\n', stdout); |
| 604 |
greg |
2.1 |
} |
