| 9 |  | */ | 
| 10 |  |  | 
| 11 |  | #include "glare.h" | 
| 12 | + | #include "linregr.h" | 
| 13 |  |  | 
| 14 |  | #define vcont(vd,vu)    ((vu)-(vd)<=SEPS) | 
| 15 |  | #define hcont(s1,s2)    ((s1)->r-(s2)->l>=-SEPS&&(s2)->r-(s1)->l>=-SEPS) | 
| 50 |  | memerr("view span brightness buffer"); | 
| 51 |  | for (v = vsize; v >= -vsize; v--) { | 
| 52 |  | close_sources(v); | 
| 53 | + | #ifndef DEBUG | 
| 54 | + | if (verbose) | 
| 55 | + | fprintf(stderr, "%s: analyzing... %3ld%%\r", | 
| 56 | + | progname, 100L*(vsize-v)/(2*vsize)); | 
| 57 | + | #endif | 
| 58 |  | getviewspan(v, spanbr); | 
| 59 |  | left = hsize + 1; | 
| 60 |  | for (h = -hsize; h <= hsize; h++) { | 
| 73 |  | addsrcspan(newspan(left,h,v,spanbr)); | 
| 74 |  | left = hsize + 1; | 
| 75 |  | } | 
| 76 | < | addindirect(h, spanbr[h+hsize]); | 
| 76 | > | addindirect(h, v, spanbr[h+hsize]); | 
| 77 |  | } | 
| 78 |  | } | 
| 79 |  | if (left < h) | 
| 80 |  | addsrcspan(newspan(left,h,v,spanbr)); | 
| 81 |  | } | 
| 76 | – | close_allsrcs(); | 
| 82 |  | free((char *)spanbr); | 
| 83 | + | close_allsrcs(); | 
| 84 |  | } | 
| 85 |  |  | 
| 86 |  |  | 
| 87 | < | addindirect(h, br)              /* add brightness to indirect illuminances */ | 
| 88 | < | int     h; | 
| 87 | > | addindirect(h, v, br)           /* add brightness to indirect illuminances */ | 
| 88 | > | int     h, v; | 
| 89 |  | double  br; | 
| 90 |  | { | 
| 91 |  | double  tanb, d; | 
| 92 | + | int     hl; | 
| 93 |  | register int    i; | 
| 94 |  |  | 
| 95 | < | if (h <= -hlim) {                       /* left region */ | 
| 96 | < | d = (double)(h+hlim)/sampdens; | 
| 97 | < | if (d <= -1.0+FTINY) | 
| 95 | > | hl = hlim(v); | 
| 96 | > | if (h <= -hl) {                 /* left region */ | 
| 97 | > | d = (double)(-h-hl)/sampdens; | 
| 98 | > | if (d >= 1.0-FTINY) | 
| 99 |  | return; | 
| 100 |  | tanb = d/sqrt(1.0-d*d); | 
| 101 |  | for (i = 0; i < nglardirs; i++) { | 
| 102 |  | d = indirect[i].lcos - tanb*indirect[i].lsin; | 
| 103 |  | if (d > 0.0) { | 
| 104 |  | indirect[i].sum += d * br; | 
| 105 | < | indirect[i].n++; | 
| 105 | > | indirect[i].n += d; | 
| 106 |  | } | 
| 107 |  | } | 
| 108 |  | return; | 
| 109 |  | } | 
| 110 | < | if (h >= hlim) {                        /* right region */ | 
| 111 | < | d = (double)(h-hlim)/sampdens; | 
| 112 | < | if (d >= 1.0-FTINY) | 
| 110 | > | if (h >= hl) {                  /* right region */ | 
| 111 | > | d = (double)(-h+hl)/sampdens; | 
| 112 | > | if (d <= -1.0+FTINY) | 
| 113 |  | return; | 
| 114 |  | tanb = d/sqrt(1.0-d*d); | 
| 115 |  | for (i = 0; i < nglardirs; i++) { | 
| 116 |  | d = indirect[i].rcos - tanb*indirect[i].rsin; | 
| 117 |  | if (d > 0.0) { | 
| 118 |  | indirect[i].sum += d * br; | 
| 119 | < | indirect[i].n++; | 
| 119 | > | indirect[i].n += d; | 
| 120 |  | } | 
| 121 |  | } | 
| 122 |  | return; | 
| 123 |  | } | 
| 124 |  | /* central region */ | 
| 125 |  | for (i = 0; i < nglardirs; i++) { | 
| 126 | < | d = cos(h_theta(h) - indirect[i].theta); | 
| 126 | > | d = cos(h_theta(h,v) - indirect[i].theta); | 
| 127 |  | if (d > 0.0) { | 
| 128 |  | indirect[i].sum += d * br; | 
| 129 | < | indirect[i].n++; | 
| 129 | > | indirect[i].n += d; | 
| 130 |  | } | 
| 131 |  | } | 
| 132 |  | } | 
| 143 |  | progname); | 
| 144 |  | brsum = 0.0; | 
| 145 |  | nsamps = 0; | 
| 146 | < | for (v = vsize; v >= -vsize; v -= TSAMPSTEP) | 
| 146 | > | for (v = vsize; v >= -vsize; v -= TSAMPSTEP) { | 
| 147 |  | for (h = -hsize; h <= hsize; h += TSAMPSTEP) { | 
| 148 |  | if ((br = getviewpix(h, v)) < 0.0) | 
| 149 |  | continue; | 
| 150 |  | brsum += br; | 
| 151 |  | nsamps++; | 
| 152 |  | } | 
| 153 | + | } | 
| 154 |  | if (nsamps == 0) { | 
| 155 |  | fprintf(stderr, "%s: no viewable scene!\n", progname); | 
| 156 |  | exit(1); | 
| 161 |  | exit(1); | 
| 162 |  | } | 
| 163 |  | if (verbose) { | 
| 164 | + | #ifdef DEBUG | 
| 165 |  | pict_stats(); | 
| 166 | + | #endif | 
| 167 |  | fprintf(stderr, | 
| 168 |  | "%s: threshold set to %f cd/m2 from %d samples\n", | 
| 169 |  | progname, threshold, nsamps); | 
| 186 |  | if (cs == NULL) | 
| 187 |  | cs = this; | 
| 188 |  | else { | 
| 178 | – | mergesource(cs, this); | 
| 189 |  | last->next = this->next; | 
| 190 | < | free((char *)this); | 
| 190 | > | mergesource(cs, this); | 
| 191 |  | this = last; | 
| 192 |  | } | 
| 193 |  | break; | 
| 199 |  | cs = (struct source *)malloc(sizeof(struct source)); | 
| 200 |  | if (cs == NULL) | 
| 201 |  | memerr("source records"); | 
| 202 | + | cs->dom = 0.0; | 
| 203 |  | cs->first = NULL; | 
| 204 |  | cs->next = curlist; | 
| 205 |  | curlist = cs; | 
| 230 |  | if (prev->next == NULL) | 
| 231 |  | prev->next = alp; | 
| 232 |  | sp->first = head.next; | 
| 233 | < | ap->first = NULL; | 
| 233 | > | if (ap->dom > 0.0 && sp->dom > 0.0) {           /* sources are done */ | 
| 234 | > | sp->dir[0] *= sp->dom; | 
| 235 | > | sp->dir[1] *= sp->dom; | 
| 236 | > | sp->dir[2] *= sp->dom; | 
| 237 | > | fvsum(sp->dir, sp->dir, ap->dir, ap->dom); | 
| 238 | > | normalize(sp->dir); | 
| 239 | > | sp->brt = (sp->brt*sp->dom + ap->brt*ap->dom) | 
| 240 | > | / (sp->dom + ap->dom); | 
| 241 | > | } | 
| 242 | > | free((char *)ap); | 
| 243 |  | } | 
| 244 |  |  | 
| 245 |  |  | 
| 276 |  | } | 
| 277 |  |  | 
| 278 |  |  | 
| 279 | + | struct srcspan * | 
| 280 | + | splitspan(sso, h, v, m)         /* divide source span at point */ | 
| 281 | + | register struct srcspan *sso; | 
| 282 | + | double  h, v, m; | 
| 283 | + | { | 
| 284 | + | register struct srcspan *ssn; | 
| 285 | + | double  d; | 
| 286 | + | int     hs; | 
| 287 | + |  | 
| 288 | + | d = h - m*(sso->v - v); | 
| 289 | + | hs = d < 0. ? d-.5 : d+.5; | 
| 290 | + | if (sso->l >= hs) | 
| 291 | + | return(NULL); | 
| 292 | + | if (sso->r <= hs) | 
| 293 | + | return(sso); | 
| 294 | + | /* need to split it */ | 
| 295 | + | ssn = (struct srcspan *)malloc(sizeof(struct srcspan)); | 
| 296 | + | if (ssn == NULL) | 
| 297 | + | memerr("source spans in splitspan"); | 
| 298 | + | ssn->brsum = (double)(hs - sso->l)/(sso->r - sso->l) * sso->brsum; | 
| 299 | + | sso->brsum -= ssn->brsum; | 
| 300 | + | ssn->v = sso->v; | 
| 301 | + | ssn->l = sso->l; | 
| 302 | + | ssn->r = sso->l = hs; | 
| 303 | + | return(ssn); | 
| 304 | + | } | 
| 305 | + |  | 
| 306 | + |  | 
| 307 | + | struct source * | 
| 308 | + | splitsource(so)                 /* divide source in two if it's big and long */ | 
| 309 | + | struct source   *so; | 
| 310 | + | { | 
| 311 | + | LRSUM   lr; | 
| 312 | + | LRLIN   fit; | 
| 313 | + | register struct srcspan *ss, *ssn; | 
| 314 | + | struct srcspan  *ssl, *ssnl, head; | 
| 315 | + | int     h; | 
| 316 | + | double  mh, mv; | 
| 317 | + | struct source   *sn; | 
| 318 | + |  | 
| 319 | + | lrclear(&lr); | 
| 320 | + | for (ss = so->first; ss != NULL; ss = ss->next) | 
| 321 | + | for (h = ss->l; h < ss->r; h++) | 
| 322 | + | lrpoint(h, ss->v, &lr); | 
| 323 | + | if ((double)lr.n/(sampdens*sampdens) < SABIG) | 
| 324 | + | return(NULL);                   /* too small */ | 
| 325 | + | if (lrfit(&fit, &lr) < 0) | 
| 326 | + | return(NULL);                   /* can't fit a line */ | 
| 327 | + | if (fit.correlation < LCORR && fit.correlation > -LCORR) | 
| 328 | + | return(NULL); | 
| 329 | + | if (verbose) | 
| 330 | + | fprintf(stderr, "%s: splitting large source\n", progname); | 
| 331 | + | mh = lrxavg(&lr); | 
| 332 | + | mv = lryavg(&lr); | 
| 333 | + | sn = (struct source *)malloc(sizeof(struct source)); | 
| 334 | + | if (sn == NULL) | 
| 335 | + | memerr("source records in splitsource"); | 
| 336 | + | sn->dom = 0.0; | 
| 337 | + | sn->first = NULL; | 
| 338 | + | ssnl = NULL; | 
| 339 | + | head.next = so->first; | 
| 340 | + | ssl = &head; | 
| 341 | + | for (ss = so->first; ss != NULL; ssl = ss, ss = ss->next) | 
| 342 | + | if ((ssn = splitspan(ss, mh, mv, fit.slope)) != NULL) { | 
| 343 | + | if (ssn == ss) {        /* remove from old */ | 
| 344 | + | ssl->next = ss->next; | 
| 345 | + | ss = ssl; | 
| 346 | + | } | 
| 347 | + | if (ssnl == NULL)       /* add to new */ | 
| 348 | + | sn->first = ssn; | 
| 349 | + | else | 
| 350 | + | ssnl->next = ssn; | 
| 351 | + | ssn->next = NULL; | 
| 352 | + | ssnl = ssn; | 
| 353 | + | } | 
| 354 | + | so->first = head.next; | 
| 355 | + | return(sn); | 
| 356 | + | } | 
| 357 | + |  | 
| 358 | + |  | 
| 359 |  | donesource(sp)                  /* finished with this source */ | 
| 360 |  | register struct source  *sp; | 
| 361 |  | { | 
| 362 | < | FVECT   dthis, dright; | 
| 362 | > | struct source   *newsrc; | 
| 363 |  | register struct srcspan *ss; | 
| 364 |  | int     h, n; | 
| 365 | < | double  d; | 
| 365 | > | double  hsum, vsum, d; | 
| 366 |  |  | 
| 367 | + | while ((newsrc = splitsource(sp)) != NULL)      /* split it? */ | 
| 368 | + | donesource(newsrc); | 
| 369 |  | sp->dom = 0.0; | 
| 370 | < | sp->dir[0] = sp->dir[1] = sp->dir[2] = 0.0; | 
| 370 | > | hsum = vsum = 0.0; | 
| 371 |  | sp->brt = 0.0; | 
| 372 |  | n = 0; | 
| 373 |  | for (ss = sp->first; ss != NULL; ss = ss->next) { | 
| 374 |  | sp->brt += ss->brsum; | 
| 375 |  | n += ss->r - ss->l; | 
| 376 | < | if (compdir(dright, ss->r, ss->v) < 0) | 
| 377 | < | compdir(dright, ss->r-2, ss->v); | 
| 378 | < | for (h = ss->r-1; h >= ss->l; h--) | 
| 379 | < | if (compdir(dthis, h, ss->v) == 0) { | 
| 380 | < | d = dist2(dthis, dright); | 
| 381 | < | fvsum(sp->dir, sp->dir, dthis, d); | 
| 280 | < | sp->dom += d; | 
| 281 | < | VCOPY(dright, dthis); | 
| 282 | < | } | 
| 283 | < | free((char *)ss); | 
| 376 | > | for (h = ss->l; h < ss->r; h++) { | 
| 377 | > | d = pixsize(h, ss->v); | 
| 378 | > | hsum += d*h; | 
| 379 | > | vsum += d*ss->v; | 
| 380 | > | sp->dom += d; | 
| 381 | > | } | 
| 382 |  | } | 
| 285 | – | sp->first = NULL; | 
| 383 |  | sp->brt /= (double)n; | 
| 384 | < | normalize(sp->dir); | 
| 384 | > | compdir(sp->dir, (int)(hsum/sp->dom), (int)(vsum/sp->dom)); | 
| 385 | > | freespans(sp); | 
| 386 |  | sp->next = donelist; | 
| 387 |  | donelist = sp; | 
| 388 |  | if (verbose) | 
| 389 |  | fprintf(stderr, | 
| 390 | < | "%s: found source at (%.3f,%.3f,%.3f), dw %.5f, br %.1f (%d samps)\n", | 
| 390 | > | "%s: source at [%.3f,%.3f,%.3f], dw %.5f, br %.1f (%d samps)\n", | 
| 391 |  | progname, sp->dir[0], sp->dir[1], sp->dir[2], | 
| 392 |  | sp->dom, sp->brt, n); | 
| 393 | + | } | 
| 394 | + |  | 
| 395 | + |  | 
| 396 | + | struct source * | 
| 397 | + | findbuddy(s, l)                 /* find close enough source to s in l*/ | 
| 398 | + | register struct source  *s, *l; | 
| 399 | + | { | 
| 400 | + | struct source   *bestbuddy = NULL; | 
| 401 | + | double  d, r, mindist = MAXBUDDY; | 
| 402 | + |  | 
| 403 | + | r = sqrt(s->dom/PI); | 
| 404 | + | for ( ; l != NULL; l = l->next) { | 
| 405 | + | d = sqrt(dist2(l->dir, s->dir)) - sqrt(l->dom/PI) - r; | 
| 406 | + | if (d < mindist) { | 
| 407 | + | bestbuddy = l; | 
| 408 | + | mindist = d; | 
| 409 | + | } | 
| 410 | + | } | 
| 411 | + | return(bestbuddy); | 
| 412 | + | } | 
| 413 | + |  | 
| 414 | + |  | 
| 415 | + | absorb_specks()                 /* eliminate too-small sources */ | 
| 416 | + | { | 
| 417 | + | struct source   head, *buddy; | 
| 418 | + | register struct source  *last, *this; | 
| 419 | + |  | 
| 420 | + | if (verbose) | 
| 421 | + | fprintf(stderr, "%s: absorbing small sources...\n", progname); | 
| 422 | + | head.next = donelist; | 
| 423 | + | last = &head; | 
| 424 | + | for (this = head.next; this != NULL; this = this->next) | 
| 425 | + | if (TOOSMALL(this)) { | 
| 426 | + | last->next = this->next; | 
| 427 | + | buddy = findbuddy(this, head.next); | 
| 428 | + | if (buddy != NULL) | 
| 429 | + | mergesource(buddy, this); | 
| 430 | + | else | 
| 431 | + | absorb(this); | 
| 432 | + | this = last; | 
| 433 | + | } else | 
| 434 | + | last = this; | 
| 435 | + | donelist = head.next; | 
| 436 | + | } | 
| 437 | + |  | 
| 438 | + |  | 
| 439 | + | absorb(s)                       /* absorb a source into indirect */ | 
| 440 | + | register struct source  *s; | 
| 441 | + | { | 
| 442 | + | FVECT   dir; | 
| 443 | + | double  d; | 
| 444 | + | register int    i; | 
| 445 | + |  | 
| 446 | + | for (i = 0; i < nglardirs; i++) { | 
| 447 | + | spinvector(dir, ourview.vdir, ourview.vup, indirect[i].theta); | 
| 448 | + | d = DOT(dir,s->dir)*s->dom*(sampdens*sampdens); | 
| 449 | + | if (d <= 0.0) | 
| 450 | + | continue; | 
| 451 | + | indirect[i].sum += d * s->brt; | 
| 452 | + | indirect[i].n += d; | 
| 453 | + | } | 
| 454 | + | freespans(s); | 
| 455 | + | free((char *)s); | 
| 456 | + | } | 
| 457 | + |  | 
| 458 | + |  | 
| 459 | + | freespans(sp)                   /* free spans associated with source */ | 
| 460 | + | struct source   *sp; | 
| 461 | + | { | 
| 462 | + | register struct srcspan *ss; | 
| 463 | + |  | 
| 464 | + | for (ss = sp->first; ss != NULL; ss = ss->next) | 
| 465 | + | free((char *)ss); | 
| 466 | + | sp->first = NULL; | 
| 467 |  | } |