| 12 |  |  | 
| 13 |  | #include "glare.h" | 
| 14 |  |  | 
| 15 | < | #define FEQ(a,b)        ((a)-(b)<=FTINY&&(a)-(b)<=FTINY) | 
| 15 | > | #define FEQ(a,b)        ((a)-(b)<=FTINY&&(b)-(a)<=FTINY) | 
| 16 |  | #define VEQ(v1,v2)      (FEQ((v1)[0],(v2)[0])&&FEQ((v1)[1],(v2)[1]) \ | 
| 17 |  | &&FEQ((v1)[2],(v2)[2])) | 
| 18 |  |  | 
| 19 | < | char    *rtargv[32] = {"rtrace", "-h", "-ov", "-fff"}; | 
| 19 | > | char    *rtargv[32] = {"rtrace", "-h-", "-ov", "-fff"}; | 
| 20 |  | int     rtargc = 4; | 
| 21 |  |  | 
| 22 | < | VIEW    ourview;                        /* our view */ | 
| 22 | > | VIEW    ourview = STDVIEW;              /* our view */ | 
| 23 |  | VIEW    pictview = STDVIEW;             /* picture view */ | 
| 24 |  | VIEW    leftview, rightview;            /* leftmost and rightmost views */ | 
| 25 |  |  | 
| 29 |  | int     verbose = 0;                    /* verbose reporting */ | 
| 30 |  | char    *progname;                      /* global argv[0] */ | 
| 31 |  |  | 
| 32 | + | double  threshold = 0.;                 /* glare threshold */ | 
| 33 | + |  | 
| 34 | + | int     sampdens = SAMPDENS;            /* sample density */ | 
| 35 |  | ANGLE   glarang[180] = {AEND};          /* glare calculation angles */ | 
| 36 |  | int     nglarangs = 0; | 
| 37 |  | double  maxtheta;                       /* maximum angle (in radians) */ | 
| 38 |  | int     hsize;                          /* horizontal size */ | 
| 36 | – | int     hlim;                           /* central limit of horizontal */ | 
| 39 |  |  | 
| 40 |  | struct illum    *indirect;              /* array of indirect illuminances */ | 
| 41 |  |  | 
| 42 | + | long    npixinvw;                       /* number of pixels in view */ | 
| 43 | + | long    npixmiss;                       /* number of pixels missed */ | 
| 44 |  |  | 
| 45 | + |  | 
| 46 |  | main(argc, argv) | 
| 47 |  | int     argc; | 
| 48 |  | char    *argv[]; | 
| 49 |  | { | 
| 50 | + | int     combine = 1; | 
| 51 |  | int     gotview = 0; | 
| 52 |  | int     rval, i; | 
| 53 |  | char    *err; | 
| 62 |  | continue; | 
| 63 |  | } | 
| 64 |  | switch (argv[i][1]) { | 
| 65 | + | case 't': | 
| 66 | + | threshold = atof(argv[++i]); | 
| 67 | + | break; | 
| 68 | + | case 'r': | 
| 69 | + | sampdens = atoi(argv[++i])/2; | 
| 70 | + | break; | 
| 71 |  | case 'v': | 
| 72 |  | if (argv[i][2] == '\0') { | 
| 73 |  | verbose++; | 
| 75 |  | } | 
| 76 |  | if (argv[i][2] != 'f') | 
| 77 |  | goto userr; | 
| 78 | < | rval = viewfile(argv[++i], &ourview, 0, 0); | 
| 78 | > | rval = viewfile(argv[++i], &ourview, NULL); | 
| 79 |  | if (rval < 0) { | 
| 80 |  | fprintf(stderr, | 
| 81 |  | "%s: cannot open view file \"%s\"\n", | 
| 101 |  | case 'p': | 
| 102 |  | picture = argv[++i]; | 
| 103 |  | break; | 
| 104 | + | case 'c': | 
| 105 | + | combine = !combine; | 
| 106 | + | break; | 
| 107 |  | case 'd': | 
| 108 | + | if (argv[i][2] == 'v') { | 
| 109 | + | rtargv[rtargc++] = argv[i]; | 
| 110 | + | break; | 
| 111 | + | } | 
| 112 | + | /* FALL THROUGH */ | 
| 113 |  | case 'l': | 
| 114 | + | case 's': | 
| 115 | + | case 'P': | 
| 116 |  | rtargv[rtargc++] = argv[i]; | 
| 117 |  | rtargv[rtargc++] = argv[++i]; | 
| 118 |  | break; | 
| 119 | + | case 'w': | 
| 120 | + | rtargv[rtargc++] = argv[i]; | 
| 121 | + | break; | 
| 122 |  | case 'a': | 
| 123 |  | rtargv[rtargc++] = argv[i]; | 
| 124 |  | if (argv[i][2] == 'v') { | 
| 140 |  | } | 
| 141 |  | /* get view */ | 
| 142 |  | if (picture != NULL) { | 
| 143 | < | rval = viewfile(picture, &pictview, 0, 0); | 
| 143 | > | rval = viewfile(picture, &pictview, NULL); | 
| 144 |  | if (rval < 0) { | 
| 145 |  | fprintf(stderr, "%s: cannot open picture file \"%s\"\n", | 
| 146 |  | progname, picture); | 
| 189 |  | exit(1); | 
| 190 |  | } | 
| 191 |  | init();                                 /* initialize program */ | 
| 192 | < | comp_thresh();                          /* compute glare threshold */ | 
| 192 | > | if (threshold <= FTINY) | 
| 193 | > | comp_thresh();                  /* compute glare threshold */ | 
| 194 |  | analyze();                              /* analyze view */ | 
| 195 | + | if (combine) | 
| 196 | + | absorb_specks();                /* eliminate tiny sources */ | 
| 197 |  | cleanup();                              /* tidy up */ | 
| 198 |  | /* print header */ | 
| 199 |  | printargs(argc, argv, stdout); | 
| 200 |  | fputs(VIEWSTR, stdout); | 
| 201 |  | fprintview(&ourview, stdout); | 
| 202 | < | printf("\n\n"); | 
| 202 | > | printf("\n"); | 
| 203 | > | fputformat("ascii", stdout); | 
| 204 | > | printf("\n"); | 
| 205 |  | printsources();                         /* print glare sources */ | 
| 206 |  | printillum();                           /* print illuminances */ | 
| 207 |  | exit(0); | 
| 213 |  | } | 
| 214 |  |  | 
| 215 |  |  | 
| 216 | + | int | 
| 217 | + | angcmp(ap1, ap2)                /* compare two angles */ | 
| 218 | + | ANGLE   *ap1, *ap2; | 
| 219 | + | { | 
| 220 | + | register int    a1, a2; | 
| 221 | + |  | 
| 222 | + | a1 = *ap1; | 
| 223 | + | a2 = *ap2; | 
| 224 | + | if (a1 == a2) { | 
| 225 | + | fprintf(stderr, "%s: duplicate glare angle (%d)\n", | 
| 226 | + | progname, a1); | 
| 227 | + | exit(1); | 
| 228 | + | } | 
| 229 | + | return(a1-a2); | 
| 230 | + | } | 
| 231 | + |  | 
| 232 | + |  | 
| 233 |  | init()                          /* initialize global variables */ | 
| 234 |  | { | 
| 235 |  | double  d; | 
| 241 |  | /* set direction vectors */ | 
| 242 |  | for (i = 0; glarang[i] != AEND; i++) | 
| 243 |  | ; | 
| 244 | < | if (i > 0 && glarang[0] <= 0 || glarang[i-1] >= 180) { | 
| 245 | < | fprintf(stderr, "%s: glare angles must be between 1 and 179\n", | 
| 244 | > | qsort(glarang, i, sizeof(ANGLE), angcmp); | 
| 245 | > | if (i > 0 && (glarang[0] <= 0 || glarang[i-1] > 180)) { | 
| 246 | > | fprintf(stderr, "%s: glare angles must be between 1 and 180\n", | 
| 247 |  | progname); | 
| 248 |  | exit(1); | 
| 249 |  | } | 
| 250 |  | nglarangs = i; | 
| 251 |  | /* nglardirs = 2*nglarangs + 1; */ | 
| 252 | < | /* vsize = SAMPDENS; */ | 
| 252 | > | /* vsize = sampdens - 1; */ | 
| 253 |  | if (nglarangs > 0) | 
| 254 |  | maxtheta = (PI/180.)*glarang[nglarangs-1]; | 
| 255 |  | else | 
| 256 |  | maxtheta = 0.0; | 
| 257 | < | hlim = SAMPDENS*maxtheta; | 
| 258 | < | hsize = SAMPDENS + hlim; | 
| 259 | < | if (hsize > (int)(PI*SAMPDENS)) | 
| 212 | < | hsize = PI*SAMPDENS; | 
| 257 | > | hsize = hlim(0) + sampdens - 1; | 
| 258 | > | if (hsize > (int)(PI*sampdens)) | 
| 259 | > | hsize = PI*sampdens; | 
| 260 |  | indirect = (struct illum *)calloc(nglardirs, sizeof(struct illum)); | 
| 261 |  | if (indirect == NULL) | 
| 262 |  | memerr("indirect illuminances"); | 
| 263 | + | npixinvw = npixmiss = 0L; | 
| 264 |  | copystruct(&leftview, &ourview); | 
| 265 |  | copystruct(&rightview, &ourview); | 
| 266 |  | spinvector(leftview.vdir, ourview.vdir, ourview.vup, maxtheta); | 
| 269 |  | setview(&rightview); | 
| 270 |  | indirect[nglarangs].lcos = | 
| 271 |  | indirect[nglarangs].rcos = cos(maxtheta); | 
| 272 | < | indirect[nglarangs].lsin = | 
| 273 | < | -(indirect[nglarangs].rsin = sin(maxtheta)); | 
| 272 | > | indirect[nglarangs].rsin = | 
| 273 | > | -(indirect[nglarangs].lsin = sin(maxtheta)); | 
| 274 |  | indirect[nglarangs].theta = 0.0; | 
| 275 |  | for (i = 0; i < nglarangs; i++) { | 
| 276 |  | d = (glarang[nglarangs-1] - glarang[i])*(PI/180.); | 
| 277 |  | indirect[nglarangs-i-1].lcos = | 
| 278 |  | indirect[nglarangs+i+1].rcos = cos(d); | 
| 279 | < | indirect[nglarangs-i-1].lsin = | 
| 280 | < | -(indirect[nglarangs+i+1].rsin = sin(d)); | 
| 279 | > | indirect[nglarangs+i+1].rsin = | 
| 280 | > | -(indirect[nglarangs-i-1].lsin = sin(d)); | 
| 281 |  | d = (glarang[nglarangs-1] + glarang[i])*(PI/180.); | 
| 282 |  | indirect[nglarangs-i-1].rcos = | 
| 283 |  | indirect[nglarangs+i+1].lcos = cos(d); | 
| 284 | < | indirect[nglarangs+i+1].lsin = | 
| 285 | < | -(indirect[nglarangs-i-1].rsin = sin(d)); | 
| 286 | < | indirect[nglarangs-i-1].theta = -(PI/180.)*glarang[i]; | 
| 287 | < | indirect[nglarangs+i+1].theta = (PI/180.)*glarang[i]; | 
| 284 | > | indirect[nglarangs-i-1].rsin = | 
| 285 | > | -(indirect[nglarangs+i+1].lsin = sin(d)); | 
| 286 | > | indirect[nglarangs-i-1].theta = (PI/180.)*glarang[i]; | 
| 287 | > | indirect[nglarangs+i+1].theta = -(PI/180.)*glarang[i]; | 
| 288 |  | } | 
| 289 |  | /* open picture */ | 
| 290 |  | if (picture != NULL) { | 
| 309 |  | cleanup()                               /* close files, wait for children */ | 
| 310 |  | { | 
| 311 |  | if (verbose) | 
| 312 | < | fprintf(stderr, "%s: cleaning up...\n", progname); | 
| 312 | > | fprintf(stderr, "%s: cleaning up...        \n", progname); | 
| 313 |  | if (picture != NULL) | 
| 314 |  | close_pict(); | 
| 315 |  | if (octree != NULL) | 
| 316 |  | done_rtrace(); | 
| 317 | + | if (npixinvw < 100*npixmiss) | 
| 318 | + | fprintf(stderr, "%s: warning -- missing %d%% of samples\n", | 
| 319 | + | progname, (int)(100L*npixmiss/npixinvw)); | 
| 320 |  | } | 
| 321 |  |  | 
| 322 |  |  | 
| 324 |  | FVECT   vd; | 
| 325 |  | int     x, y; | 
| 326 |  | { | 
| 327 | + | int     hl; | 
| 328 |  | FVECT   org;                    /* dummy variable */ | 
| 329 |  |  | 
| 330 | < | if (x <= -hlim)                 /* left region */ | 
| 330 | > | hl = hlim(y); | 
| 331 | > | if (x <= -hl) {                 /* left region */ | 
| 332 | > | if (x <= -hl-sampdens) | 
| 333 | > | return(-1); | 
| 334 |  | return(viewray(org, vd, &leftview, | 
| 335 | < | (x+hlim)/(2.*SAMPDENS)+.5, | 
| 336 | < | y/(2.*SAMPDENS)+.5)); | 
| 337 | < | if (x >= hlim)                  /* right region */ | 
| 335 | > | (double)(x+hl)/(2*sampdens)+.5, | 
| 336 | > | (double)y/(2*sampdens)+.5)); | 
| 337 | > | } | 
| 338 | > | if (x >= hl) {                  /* right region */ | 
| 339 | > | if (x >= hl+sampdens) | 
| 340 | > | return(-1); | 
| 341 |  | return(viewray(org, vd, &rightview, | 
| 342 | < | (x-hlim)/(2.*SAMPDENS)+.5, | 
| 343 | < | y/(2.*SAMPDENS)+.5)); | 
| 344 | < | /* central region */ | 
| 345 | < | if (viewray(org, vd, &ourview, .5, y/(2.*SAMPDENS)+.5) < 0) | 
| 342 | > | (double)(x-hl)/(2*sampdens)+.5, | 
| 343 | > | (double)y/(2*sampdens)+.5)); | 
| 344 | > | } | 
| 345 | > | /* central region */ | 
| 346 | > | if (viewray(org, vd, &ourview, .5, (double)y/(2*sampdens)+.5) < 0) | 
| 347 |  | return(-1); | 
| 348 | < | spinvector(vd, vd, ourview.vup, h_theta(x)); | 
| 348 | > | spinvector(vd, vd, ourview.vup, h_theta(x,y)); | 
| 349 |  | return(0); | 
| 350 |  | } | 
| 351 |  |  | 
| 352 |  |  | 
| 353 | < | spinvector(vres, vorig, vnorm, theta)   /* rotate vector around normal */ | 
| 354 | < | FVECT  vres, vorig, vnorm; | 
| 355 | < | double  theta; | 
| 353 | > | double | 
| 354 | > | pixsize(x, y)           /* return the solid angle of pixel at (x,y) */ | 
| 355 | > | int     x, y; | 
| 356 |  | { | 
| 357 | < | extern double  sin(), cos(); | 
| 358 | < | double  sint, cost, dotp; | 
| 359 | < | FVECT  vperp; | 
| 360 | < | register int  i; | 
| 361 | < |  | 
| 362 | < | sint = sin(theta); | 
| 363 | < | cost = cos(theta); | 
| 364 | < | dotp = DOT(vorig, vnorm); | 
| 365 | < | fcross(vperp, vnorm, vorig); | 
| 366 | < | for (i = 0; i < 3; i++) | 
| 367 | < | vres[i] = vnorm[i]*dotp*(1.-cost) + | 
| 368 | < | vorig[i]*cost + vperp[i]*sint; | 
| 357 | > | register int    hl, xo; | 
| 358 | > | double  disc; | 
| 359 | > |  | 
| 360 | > | hl = hlim(y); | 
| 361 | > | if (x < -hl) | 
| 362 | > | xo = x+hl; | 
| 363 | > | else if (x > hl) | 
| 364 | > | xo = x-hl; | 
| 365 | > | else | 
| 366 | > | xo = 0; | 
| 367 | > | disc = 1. - (double)((long)xo*xo + (long)y*y)/((long)sampdens*sampdens); | 
| 368 | > | if (disc <= FTINY*FTINY) | 
| 369 | > | return(0.); | 
| 370 | > | return(1./(sampdens*sampdens*sqrt(disc))); | 
| 371 |  | } | 
| 372 |  |  | 
| 373 |  |  | 
| 374 |  | memerr(s)                       /* malloc failure */ | 
| 375 |  | char    *s; | 
| 376 |  | { | 
| 377 | < | fprintf(stderr, "%s: out of memory for %s\n", s); | 
| 377 | > | fprintf(stderr, "%s: out of memory for %s\n", progname, s); | 
| 378 |  | exit(1); | 
| 379 |  | } | 
| 380 |  |  | 
| 398 |  |  | 
| 399 |  | printf("BEGIN indirect illuminance\n"); | 
| 400 |  | for (i = 0; i < nglardirs; i++) | 
| 401 | < | printf("\t%.0f\t%f\n", (180.0/PI)*indirect[i].theta, | 
| 402 | < | PI * indirect[i].sum / (double)indirect[i].n); | 
| 401 | > | if (indirect[i].n > FTINY) | 
| 402 | > | printf("\t%.0f\t%f\n", (180.0/PI)*indirect[i].theta, | 
| 403 | > | PI * indirect[i].sum / indirect[i].n); | 
| 404 |  | printf("END indirect illuminance\n"); | 
| 405 |  | } |