| 1 | #ifndef lint | 
| 2 | static const char       RCSid[] = "$Id: findglare.c,v 2.17 2023/01/15 16:38:10 greg Exp $"; | 
| 3 | #endif | 
| 4 | /* | 
| 5 | * Find glare sources in a scene or image. | 
| 6 | * | 
| 7 | *      Greg Ward       March 1991 | 
| 8 | */ | 
| 9 |  | 
| 10 | #include "glare.h" | 
| 11 |  | 
| 12 | char    *rtargv[64] = {"rtrace", "-h-", "-ov", "-fff", "-ld-", "-i-", "-I-"}; | 
| 13 | int     rtargc = 7; | 
| 14 |  | 
| 15 | VIEW    ourview = STDVIEW;              /* our view */ | 
| 16 | VIEW    pictview = STDVIEW;             /* picture view */ | 
| 17 | VIEW    leftview, rightview;            /* leftmost and rightmost views */ | 
| 18 |  | 
| 19 | char    *picture = NULL;                /* picture file name */ | 
| 20 | char    *octree = NULL;                 /* octree file name */ | 
| 21 |  | 
| 22 | int     verbose = 0;                    /* verbose reporting */ | 
| 23 |  | 
| 24 | double  threshold = 0.;                 /* glare threshold */ | 
| 25 |  | 
| 26 | int     sampdens = SAMPDENS;            /* sample density */ | 
| 27 | ANGLE   glarang[180] = {AEND};          /* glare calculation angles */ | 
| 28 | int     nglarangs = 0; | 
| 29 | double  maxtheta;                       /* maximum angle (in radians) */ | 
| 30 | int     hsize;                          /* horizontal size */ | 
| 31 |  | 
| 32 | struct illum    *indirect;              /* array of indirect illuminances */ | 
| 33 |  | 
| 34 | long    npixinvw;                       /* number of pixels in view */ | 
| 35 | long    npixmiss;                       /* number of pixels missed */ | 
| 36 |  | 
| 37 | static int angcmp(const void *ap1, const void *ap2); | 
| 38 | static void init(void); | 
| 39 | static void cleanup(void); | 
| 40 | static void printsources(void); | 
| 41 | static void printillum(void); | 
| 42 |  | 
| 43 |  | 
| 44 |  | 
| 45 | int | 
| 46 | main( | 
| 47 | int     argc, | 
| 48 | char    *argv[] | 
| 49 | ) | 
| 50 | { | 
| 51 | int     combine = 1; | 
| 52 | int     gotview = 0; | 
| 53 | int     rval, i; | 
| 54 | char    *err; | 
| 55 | /* set global progname */ | 
| 56 | fixargv0(argv[0]); | 
| 57 | /* process options */ | 
| 58 | for (i = 1; i < argc && argv[i][0] == '-'; i++) { | 
| 59 | /* expand arguments */ | 
| 60 | while ((rval = expandarg(&argc, &argv, i)) > 0) | 
| 61 | ; | 
| 62 | if (rval < 0) { | 
| 63 | fprintf(stderr, "%s: cannot expand '%s'\n", | 
| 64 | argv[0], argv[i]); | 
| 65 | exit(1); | 
| 66 | } | 
| 67 | rval = getviewopt(&ourview, argc-i, argv+i); | 
| 68 | if (rval >= 0) { | 
| 69 | i += rval; | 
| 70 | gotview++; | 
| 71 | continue; | 
| 72 | } | 
| 73 | switch (argv[i][1]) { | 
| 74 | case 't': | 
| 75 | threshold = atof(argv[++i]); | 
| 76 | break; | 
| 77 | case 'r': | 
| 78 | sampdens = atoi(argv[++i])/2; | 
| 79 | break; | 
| 80 | case 'v': | 
| 81 | if (argv[i][2] == '\0') { | 
| 82 | verbose++; | 
| 83 | break; | 
| 84 | } | 
| 85 | if (argv[i][2] != 'f') | 
| 86 | goto userr; | 
| 87 | rval = viewfile(argv[++i], &ourview, NULL); | 
| 88 | if (rval < 0) { | 
| 89 | fprintf(stderr, | 
| 90 | "%s: cannot open view file \"%s\"\n", | 
| 91 | progname, argv[i]); | 
| 92 | exit(1); | 
| 93 | } else if (rval == 0) { | 
| 94 | fprintf(stderr, | 
| 95 | "%s: bad view file \"%s\"\n", | 
| 96 | progname, argv[i]); | 
| 97 | exit(1); | 
| 98 | } else | 
| 99 | gotview++; | 
| 100 | break; | 
| 101 | case 'g': | 
| 102 | if (argv[i][2] != 'a') | 
| 103 | goto userr; | 
| 104 | if (setscan(glarang, argv[++i]) < 0) { | 
| 105 | fprintf(stderr, "%s: bad angle spec \"%s\"\n", | 
| 106 | progname, argv[i]); | 
| 107 | exit(1); | 
| 108 | } | 
| 109 | break; | 
| 110 | case 'p': | 
| 111 | picture = argv[++i]; | 
| 112 | break; | 
| 113 | case 'c': | 
| 114 | combine = !combine; | 
| 115 | break; | 
| 116 | case 'd': | 
| 117 | rtargv[rtargc++] = argv[i]; | 
| 118 | if (argv[i][2] != 'v') | 
| 119 | rtargv[rtargc++] = argv[++i]; | 
| 120 | break; | 
| 121 | case 'l': | 
| 122 | if (argv[i][2] == 'd') | 
| 123 | break; | 
| 124 | /* FALL THROUGH */ | 
| 125 | case 's': | 
| 126 | case 'P': | 
| 127 | case 'n': | 
| 128 | rtargv[rtargc++] = argv[i]; | 
| 129 | rtargv[rtargc++] = argv[++i]; | 
| 130 | break; | 
| 131 | case 'w': | 
| 132 | case 'u': | 
| 133 | case 'b': | 
| 134 | rtargv[rtargc++] = argv[i]; | 
| 135 | break; | 
| 136 | case 'a': | 
| 137 | rtargv[rtargc++] = argv[i]; | 
| 138 | if (argv[i][2] == 'v') { | 
| 139 | rtargv[rtargc++] = argv[++i]; | 
| 140 | rtargv[rtargc++] = argv[++i]; | 
| 141 | } | 
| 142 | rtargv[rtargc++] = argv[++i]; | 
| 143 | break; | 
| 144 | case 'm': | 
| 145 | rtargv[rtargc++] = argv[i]; | 
| 146 | if ((argv[i][2] == 'e') | (argv[i][2] == 'a')) { | 
| 147 | rtargv[rtargc++] = argv[++i]; | 
| 148 | rtargv[rtargc++] = argv[++i]; | 
| 149 | } | 
| 150 | rtargv[rtargc++] = argv[++i]; | 
| 151 | break; | 
| 152 | default: | 
| 153 | goto userr; | 
| 154 | } | 
| 155 | } | 
| 156 | /* get octree */ | 
| 157 | if (i < argc-1) | 
| 158 | goto userr; | 
| 159 | if (i == argc-1) { | 
| 160 | rtargv[rtargc++] = octree = argv[i]; | 
| 161 | rtargv[rtargc] = NULL; | 
| 162 | } | 
| 163 | /* get view */ | 
| 164 | if (picture != NULL) { | 
| 165 | rval = viewfile(picture, &pictview, NULL); | 
| 166 | if (rval < 0) { | 
| 167 | fprintf(stderr, "%s: cannot open picture file \"%s\"\n", | 
| 168 | progname, picture); | 
| 169 | exit(1); | 
| 170 | } else if (rval == 0) { | 
| 171 | fprintf(stderr, | 
| 172 | "%s: cannot get view from picture \"%s\"\n", | 
| 173 | progname, picture); | 
| 174 | exit(1); | 
| 175 | } | 
| 176 | if (pictview.type == VT_PAR) { | 
| 177 | fprintf(stderr, "%s: %s: cannot use parallel view\n", | 
| 178 | progname, picture); | 
| 179 | exit(1); | 
| 180 | } | 
| 181 | if ((err = setview(&pictview)) != NULL) { | 
| 182 | fprintf(stderr, "%s: %s\n", picture, err); | 
| 183 | exit(1); | 
| 184 | } | 
| 185 | } | 
| 186 | if (!gotview) { | 
| 187 | if (picture == NULL) { | 
| 188 | fprintf(stderr, "%s: must have view or picture\n", | 
| 189 | progname); | 
| 190 | exit(1); | 
| 191 | } | 
| 192 | ourview = pictview; | 
| 193 | } else if (picture != NULL && !VABSEQ(ourview.vp, pictview.vp)) { | 
| 194 | fprintf(stderr, "%s: picture must have same viewpoint\n", | 
| 195 | progname); | 
| 196 | exit(1); | 
| 197 | } | 
| 198 | ourview.type = VT_HEM; | 
| 199 | ourview.horiz = ourview.vert = 180.0; | 
| 200 | ourview.hoff = ourview.voff = 0.0; | 
| 201 | fvsum(ourview.vdir, ourview.vdir, ourview.vup, | 
| 202 | -DOT(ourview.vdir,ourview.vup)); | 
| 203 | if ((err = setview(&ourview)) != NULL) { | 
| 204 | fprintf(stderr, "%s: %s\n", progname, err); | 
| 205 | exit(1); | 
| 206 | } | 
| 207 | if (octree == NULL && picture == NULL) { | 
| 208 | fprintf(stderr, | 
| 209 | "%s: must specify at least one of picture or octree\n", | 
| 210 | progname); | 
| 211 | exit(1); | 
| 212 | } | 
| 213 | init();                                 /* initialize program */ | 
| 214 | if (threshold <= FTINY) | 
| 215 | comp_thresh();                  /* compute glare threshold */ | 
| 216 | analyze();                              /* analyze view */ | 
| 217 | if (combine) | 
| 218 | absorb_specks();                /* eliminate tiny sources */ | 
| 219 | cleanup();                              /* tidy up */ | 
| 220 | /* print header */ | 
| 221 | newheader("RADIANCE", stdout); | 
| 222 | printargs(argc, argv, stdout); | 
| 223 | fputs(VIEWSTR, stdout); | 
| 224 | fprintview(&ourview, stdout); | 
| 225 | printf("\n"); | 
| 226 | fputformat("ascii", stdout); | 
| 227 | printf("\n"); | 
| 228 | printsources();                         /* print glare sources */ | 
| 229 | printillum();                           /* print illuminances */ | 
| 230 | exit(0); | 
| 231 | userr: | 
| 232 | fprintf(stderr, | 
| 233 | "Usage: %s [view options][-ga angles][-p picture][[rtrace options] octree]\n", | 
| 234 | progname); | 
| 235 | exit(1); | 
| 236 | } | 
| 237 |  | 
| 238 |  | 
| 239 | static int | 
| 240 | angcmp(         /* compare two angles */ | 
| 241 | const void      *ap1, | 
| 242 | const void      *ap2 | 
| 243 | ) | 
| 244 | { | 
| 245 | int     a1, a2; | 
| 246 |  | 
| 247 | a1 = *(ANGLE *)ap1; | 
| 248 | a2 = *(ANGLE *)ap2; | 
| 249 | if (a1 == a2) { | 
| 250 | fprintf(stderr, "%s: duplicate glare angle (%d)\n", | 
| 251 | progname, a1); | 
| 252 | exit(1); | 
| 253 | } | 
| 254 | return(a1-a2); | 
| 255 | } | 
| 256 |  | 
| 257 |  | 
| 258 | static void | 
| 259 | init(void)                              /* initialize global variables */ | 
| 260 | { | 
| 261 | double  d; | 
| 262 | int     i; | 
| 263 |  | 
| 264 | if (verbose) | 
| 265 | fprintf(stderr, "%s: initializing data structures...\n", | 
| 266 | progname); | 
| 267 | /* set direction vectors */ | 
| 268 | for (i = 0; glarang[i] != AEND; i++) | 
| 269 | ; | 
| 270 | qsort(glarang, i, sizeof(ANGLE), angcmp); | 
| 271 | if (i > 0 && (glarang[0] <= 0 || glarang[i-1] > 180)) { | 
| 272 | fprintf(stderr, "%s: glare angles must be between 1 and 180\n", | 
| 273 | progname); | 
| 274 | exit(1); | 
| 275 | } | 
| 276 | nglarangs = i; | 
| 277 | /* nglardirs = 2*nglarangs + 1; */ | 
| 278 | /* vsize = sampdens - 1; */ | 
| 279 | if (nglarangs > 0) | 
| 280 | maxtheta = (PI/180.)*glarang[nglarangs-1]; | 
| 281 | else | 
| 282 | maxtheta = 0.0; | 
| 283 | hsize = hlim(0) + sampdens - 1; | 
| 284 | if (hsize > (int)(PI*sampdens)) | 
| 285 | hsize = PI*sampdens; | 
| 286 | indirect = (struct illum *)calloc(nglardirs, sizeof(struct illum)); | 
| 287 | if (indirect == NULL) | 
| 288 | memerr("indirect illuminances"); | 
| 289 | npixinvw = npixmiss = 0L; | 
| 290 | leftview = ourview; | 
| 291 | rightview = ourview; | 
| 292 | spinvector(leftview.vdir, ourview.vdir, ourview.vup, maxtheta); | 
| 293 | spinvector(rightview.vdir, ourview.vdir, ourview.vup, -maxtheta); | 
| 294 | setview(&leftview); | 
| 295 | setview(&rightview); | 
| 296 | indirect[nglarangs].lcos = | 
| 297 | indirect[nglarangs].rcos = cos(maxtheta); | 
| 298 | indirect[nglarangs].rsin = | 
| 299 | -(indirect[nglarangs].lsin = sin(maxtheta)); | 
| 300 | indirect[nglarangs].theta = 0.0; | 
| 301 | for (i = 0; i < nglarangs; i++) { | 
| 302 | d = (glarang[nglarangs-1] - glarang[i])*(PI/180.); | 
| 303 | indirect[nglarangs-i-1].lcos = | 
| 304 | indirect[nglarangs+i+1].rcos = cos(d); | 
| 305 | indirect[nglarangs+i+1].rsin = | 
| 306 | -(indirect[nglarangs-i-1].lsin = sin(d)); | 
| 307 | d = (glarang[nglarangs-1] + glarang[i])*(PI/180.); | 
| 308 | indirect[nglarangs-i-1].rcos = | 
| 309 | indirect[nglarangs+i+1].lcos = cos(d); | 
| 310 | indirect[nglarangs-i-1].rsin = | 
| 311 | -(indirect[nglarangs+i+1].lsin = sin(d)); | 
| 312 | indirect[nglarangs-i-1].theta = (PI/180.)*glarang[i]; | 
| 313 | indirect[nglarangs+i+1].theta = -(PI/180.)*glarang[i]; | 
| 314 | } | 
| 315 | /* open picture */ | 
| 316 | if (picture != NULL) { | 
| 317 | if (verbose) | 
| 318 | fprintf(stderr, "%s: opening picture file \"%s\"...\n", | 
| 319 | progname, picture); | 
| 320 | open_pict(picture); | 
| 321 | } | 
| 322 | /* start rtrace */ | 
| 323 | if (octree != NULL) { | 
| 324 | if (verbose) { | 
| 325 | fprintf(stderr, | 
| 326 | "%s: starting luminance calculation...\n\t", | 
| 327 | progname); | 
| 328 | printargs(rtargc, rtargv, stderr); | 
| 329 | } | 
| 330 | fork_rtrace(rtargv); | 
| 331 | } | 
| 332 | } | 
| 333 |  | 
| 334 |  | 
| 335 | static void | 
| 336 | cleanup(void)                           /* close files, wait for children */ | 
| 337 | { | 
| 338 | if (verbose) | 
| 339 | fprintf(stderr, "%s: cleaning up...        \n", progname); | 
| 340 | if (picture != NULL) | 
| 341 | close_pict(); | 
| 342 | if (octree != NULL) | 
| 343 | done_rtrace(); | 
| 344 | if (npixinvw < 100*npixmiss) | 
| 345 | fprintf(stderr, "%s: warning -- missing %d%% of samples\n", | 
| 346 | progname, (int)(100L*npixmiss/npixinvw)); | 
| 347 | } | 
| 348 |  | 
| 349 |  | 
| 350 | int | 
| 351 | compdir(                        /* compute direction for x,y */ | 
| 352 | FVECT   vd, | 
| 353 | int     x, | 
| 354 | int     y | 
| 355 | ) | 
| 356 | { | 
| 357 | int     hl; | 
| 358 | FVECT   org;                    /* dummy variable */ | 
| 359 |  | 
| 360 | hl = hlim(y); | 
| 361 | if (x <= -hl) {                 /* left region */ | 
| 362 | if (x <= -hl-sampdens) | 
| 363 | return(-1); | 
| 364 | return(viewray(org, vd, &leftview, | 
| 365 | (double)(x+hl)/(2*sampdens)+.5, | 
| 366 | (double)y/(2*sampdens)+.5)); | 
| 367 | } | 
| 368 | if (x >= hl) {                  /* right region */ | 
| 369 | if (x >= hl+sampdens) | 
| 370 | return(-1); | 
| 371 | return(viewray(org, vd, &rightview, | 
| 372 | (double)(x-hl)/(2*sampdens)+.5, | 
| 373 | (double)y/(2*sampdens)+.5)); | 
| 374 | } | 
| 375 | /* central region */ | 
| 376 | if (viewray(org, vd, &ourview, .5, (double)y/(2*sampdens)+.5) < 0) | 
| 377 | return(-1); | 
| 378 | spinvector(vd, vd, ourview.vup, h_theta(x,y)); | 
| 379 | return(0); | 
| 380 | } | 
| 381 |  | 
| 382 |  | 
| 383 | double | 
| 384 | pixsize(                /* return the solid angle of pixel at (x,y) */ | 
| 385 | int     x, | 
| 386 | int     y | 
| 387 | ) | 
| 388 | { | 
| 389 | int     hl, xo; | 
| 390 | double  disc; | 
| 391 |  | 
| 392 | hl = hlim(y); | 
| 393 | if (x < -hl) | 
| 394 | xo = x+hl; | 
| 395 | else if (x > hl) | 
| 396 | xo = x-hl; | 
| 397 | else | 
| 398 | xo = 0; | 
| 399 | disc = 1. - (double)((long)xo*xo + (long)y*y)/((long)sampdens*sampdens); | 
| 400 | if (disc <= FTINY*FTINY) | 
| 401 | return(0.); | 
| 402 | return(1./(sampdens*sampdens*sqrt(disc))); | 
| 403 | } | 
| 404 |  | 
| 405 |  | 
| 406 | void | 
| 407 | memerr(                 /* malloc failure */ | 
| 408 | char    *s | 
| 409 | ) | 
| 410 | { | 
| 411 | fprintf(stderr, "%s: out of memory for %s\n", progname, s); | 
| 412 | exit(1); | 
| 413 | } | 
| 414 |  | 
| 415 |  | 
| 416 | static void | 
| 417 | printsources(void)                      /* print out glare sources */ | 
| 418 | { | 
| 419 | struct source   *sp; | 
| 420 |  | 
| 421 | printf("BEGIN glare source\n"); | 
| 422 | for (sp = donelist; sp != NULL; sp = sp->next) | 
| 423 | printf("\t%f %f %f\t%e\t%e\n", | 
| 424 | sp->dir[0], sp->dir[1], sp->dir[2], | 
| 425 | sp->dom, sp->brt); | 
| 426 | printf("END glare source\n"); | 
| 427 | } | 
| 428 |  | 
| 429 |  | 
| 430 | static void | 
| 431 | printillum(void)                        /* print out indirect illuminances */ | 
| 432 | { | 
| 433 | int     i; | 
| 434 |  | 
| 435 | printf("BEGIN indirect illuminance\n"); | 
| 436 | for (i = 0; i < nglardirs; i++) | 
| 437 | if (indirect[i].n > FTINY) | 
| 438 | printf("\t%.0f\t%f\n", (180.0/PI)*indirect[i].theta, | 
| 439 | PI * indirect[i].sum / indirect[i].n); | 
| 440 | printf("END indirect illuminance\n"); | 
| 441 | } |