Back to index

radiance  4R0+20100331
glaresrc.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: glaresrc.c,v 2.6 2006/05/29 16:47:54 greg Exp $";
00003 #endif
00004 /*
00005  * Gather samples and compute glare sources.
00006  */
00007 
00008 #include "glare.h"
00009 #include "linregr.h"
00010 
00011 #define vcont(vd,vu) ((vu)-(vd)<=SEPS)
00012 #define hcont(s1,s2) ((s1)->r-(s2)->l>=-SEPS&&(s2)->r-(s1)->l>=-SEPS)
00013 
00014 struct source *curlist = NULL;     /* current source list */
00015 struct source *donelist = NULL;    /* finished sources */
00016 
00017 void   pict_stats(void);
00018 
00019 static struct srcspan * newspan(int       l, int r, int v, float      *sb);
00020 static struct srcspan * newspan(int       l, int r, int v, float      *sb);
00021 static void addindirect(int h, int v, double     br);
00022 static void addsrcspan(struct srcspan     *nss);
00023 static void mergesource(struct source     *sp, struct source   *ap);
00024 static void close_sources(int      v);
00025 static void close_allsrcs(void);
00026 static struct srcspan * splitspan(register struct srcspan      *sso, double  h, double     v, double     m);
00027 static struct source * splitsource(struct source *so);
00028 static void donesource(register struct source    *sp);
00029 static struct source * findbuddy(register struct source *s, register struct source  *l);
00030 static void absorb(register struct source *s);
00031 static void freespans(struct source       *sp);
00032 
00033 
00034 static struct srcspan *
00035 newspan(             /* allocate a new source span */
00036        int    l,
00037        int    r,
00038        int    v,
00039        float  *sb
00040 )
00041 {
00042        register struct srcspan     *ss;
00043        register int  i;
00044 
00045        ss = (struct srcspan *)malloc(sizeof(struct srcspan));
00046        if (ss == NULL)
00047               memerr("source spans");
00048        ss->l = l;
00049        ss->r = r;
00050        ss->v = v;
00051        ss->brsum = 0.0;
00052        for (i = l; i < r; i++)
00053               ss->brsum += sb[i+hsize];
00054        return(ss);
00055 }
00056 
00057 
00058 extern void
00059 analyze(void)               /* analyze our scene */
00060 {
00061        int    h, v;
00062        int    left;
00063        float  *spanbr;
00064 
00065        spanbr = (float *)malloc((2*hsize+1)*sizeof(float));
00066        if (spanbr == NULL)
00067               memerr("view span brightness buffer");
00068        for (v = vsize; v >= -vsize; v--) {
00069               close_sources(v);
00070 #ifndef DEBUG
00071               if (verbose) {
00072                      fprintf(stderr, "%s: analyzing... %3ld%%\r",
00073                             progname, 100L*(vsize-v)/(2*vsize));
00074                      fflush(stderr);
00075               }
00076 #endif
00077               getviewspan(v, spanbr);
00078               left = hsize + 1;
00079               for (h = -hsize; h <= hsize; h++) {
00080                      if (spanbr[h+hsize] < 0.0) {       /* off view */
00081                             if (left < h) {
00082                                    addsrcspan(newspan(left,h,v,spanbr));
00083                                    left = hsize + 1;
00084                             }
00085                             continue;
00086                      }
00087                      if (spanbr[h+hsize] > threshold) { /* in source */
00088                             if (left > h)
00089                                    left = h;
00090                      } else {                    /* out of source */
00091                             if (left < h) {
00092                                    addsrcspan(newspan(left,h,v,spanbr));
00093                                    left = hsize + 1;
00094                             }
00095                             addindirect(h, v, spanbr[h+hsize]);
00096                      }
00097               }
00098               if (left < h)
00099                      addsrcspan(newspan(left,h,v,spanbr));
00100        }
00101        free((void *)spanbr);
00102        close_allsrcs();
00103 }
00104 
00105 
00106 static void
00107 addindirect(         /* add brightness to indirect illuminances */
00108        int    h,
00109        int    v,
00110        double br
00111 )
00112 {
00113        double tanb, d;
00114        int    hl;
00115        register int  i;
00116 
00117        hl = hlim(v);
00118        if (h <= -hl) {                    /* left region */
00119               d = (double)(-h-hl)/sampdens;
00120               if (d >= 1.0-FTINY)
00121                      return;
00122               tanb = d/sqrt(1.0-d*d);
00123               for (i = 0; i < nglardirs; i++) {
00124                      d = indirect[i].lcos - tanb*indirect[i].lsin;
00125                      if (d > 0.0) {
00126                             indirect[i].sum += d * br;
00127                             indirect[i].n += d;
00128                      }
00129               }
00130               return;
00131        }
00132        if (h >= hl) {                     /* right region */
00133               d = (double)(-h+hl)/sampdens;
00134               if (d <= -1.0+FTINY)
00135                      return;
00136               tanb = d/sqrt(1.0-d*d);
00137               for (i = 0; i < nglardirs; i++) {
00138                      d = indirect[i].rcos - tanb*indirect[i].rsin;
00139                      if (d > 0.0) {
00140                             indirect[i].sum += d * br;
00141                             indirect[i].n += d;
00142                      }
00143               }
00144               return;
00145        }
00146                                    /* central region */
00147        for (i = 0; i < nglardirs; i++) {
00148               d = cos(h_theta(h,v) - indirect[i].theta);
00149               if (d > 0.0) {
00150                      indirect[i].sum += d * br;
00151                      indirect[i].n += d;
00152               }
00153        }
00154 }
00155 
00156 
00157 extern void
00158 comp_thresh(void)                  /* compute glare threshold */
00159 {
00160        int    h, v;
00161        int    nsamps;
00162        double brsum, br;
00163 
00164        if (verbose)
00165               fprintf(stderr, "%s: computing glare threshold...\n",
00166                             progname);
00167        brsum = 0.0;
00168        nsamps = 0;
00169        for (v = vsize; v >= -vsize; v -= TSAMPSTEP) {
00170               for (h = -hsize; h <= hsize; h += TSAMPSTEP) {
00171                      if ((br = getviewpix(h, v)) < 0.0)
00172                             continue;
00173                      brsum += br;
00174                      nsamps++;
00175               }
00176        }
00177        if (nsamps == 0) {
00178               fprintf(stderr, "%s: no viewable scene!\n", progname);
00179               exit(1);
00180        }
00181        threshold = GLAREBR * brsum / nsamps;
00182        if (threshold <= FTINY) {
00183               fprintf(stderr, "%s: threshold zero!\n", progname);
00184               exit(1);
00185        }
00186        if (verbose) {
00187 #ifdef DEBUG
00188               pict_stats();
00189 #endif
00190               fprintf(stderr,
00191                      "%s: threshold set to %f cd/m2 from %d samples\n",
00192                             progname, threshold, nsamps);
00193        }
00194 }
00195 
00196 
00197 static void
00198 addsrcspan(                 /* add new source span to our list */
00199        struct srcspan       *nss
00200 )
00201 {
00202        struct source *last, *cs, *this;
00203        register struct srcspan     *ss;
00204 
00205        cs = NULL;
00206        for (this = curlist; this != NULL; this = this->next) {
00207               for (ss = this->first; ss != NULL; ss = ss->next) {
00208                      if (!vcont(nss->v, ss->v))
00209                             break;
00210                      if (hcont(ss, nss)) {
00211                             if (cs == NULL)
00212                                    cs = this;
00213                             else {
00214                                    last->next = this->next;
00215                                    mergesource(cs, this);
00216                                    this = last;
00217                             }
00218                             break;
00219                      }
00220               }
00221               last = this;
00222        }
00223        if (cs == NULL) {
00224               cs = (struct source *)malloc(sizeof(struct source));
00225               if (cs == NULL)
00226                      memerr("source records");
00227               cs->dom = 0.0;
00228               cs->first = NULL;
00229               cs->next = curlist;
00230               curlist = cs;
00231        }
00232        nss->next = cs->first;
00233        cs->first = nss;
00234 }
00235 
00236 
00237 static void
00238 mergesource(         /* merge source ap into source sp */
00239        struct source *sp,
00240        struct source *ap
00241 )
00242 {
00243        struct srcspan       head;
00244        register struct srcspan     *alp, *prev, *tp;
00245 
00246        head.next = sp->first;
00247        prev = &head;
00248        alp = ap->first;
00249        while (alp != NULL && prev->next != NULL) {
00250               if (prev->next->v > alp->v) {
00251                      tp = alp->next;
00252                      alp->next = prev->next;
00253                      prev->next = alp;
00254                      alp = tp;
00255               }
00256               prev = prev->next;
00257        }
00258        if (prev->next == NULL)
00259               prev->next = alp;
00260        sp->first = head.next;
00261        if (ap->dom > 0.0 && sp->dom > 0.0) {            /* sources are done */
00262               sp->dir[0] *= sp->dom;
00263               sp->dir[1] *= sp->dom;
00264               sp->dir[2] *= sp->dom;
00265               fvsum(sp->dir, sp->dir, ap->dir, ap->dom);
00266               normalize(sp->dir);
00267               sp->brt = (sp->brt*sp->dom + ap->brt*ap->dom)
00268                             / (sp->dom + ap->dom);
00269        }
00270        free((void *)ap);
00271 }
00272 
00273 
00274 static void
00275 close_sources(              /* close sources above v */
00276        int    v
00277 )
00278 {
00279        struct source head;
00280        register struct source      *last, *this;
00281 
00282        head.next = curlist;
00283        last = &head;
00284        for (this = curlist; this != NULL; this = this->next)
00285               if (!vcont(v, this->first->v)) {
00286                      last->next = this->next;
00287                      donesource(this);
00288                      this = last;
00289               } else
00290                      last = this;
00291        curlist = head.next;
00292 }
00293 
00294 
00295 static void
00296 close_allsrcs(void)                /* done with everything */
00297 {
00298        register struct source      *this, *next;
00299 
00300        this = curlist;
00301        while (this != NULL) {
00302               next = this->next;
00303               donesource(this);
00304               this = next;
00305        }
00306        curlist = NULL;
00307 }
00308 
00309 
00310 static struct srcspan *
00311 splitspan(           /* divide source span at point */
00312        register struct srcspan     *sso,
00313        double h,
00314        double v,
00315        double m
00316 )
00317 {
00318        register struct srcspan     *ssn;
00319        double d;
00320        int    hs;
00321 
00322        d = h - m*(sso->v - v);
00323        hs = d < 0. ? d-.5 : d+.5;
00324        if (sso->l >= hs)
00325               return(NULL);
00326        if (sso->r <= hs)
00327               return(sso);
00328                             /* need to split it */
00329        ssn = (struct srcspan *)malloc(sizeof(struct srcspan));
00330        if (ssn == NULL)
00331               memerr("source spans in splitspan");
00332        ssn->brsum = (double)(hs - sso->l)/(sso->r - sso->l) * sso->brsum;
00333        sso->brsum -= ssn->brsum;
00334        ssn->v = sso->v;
00335        ssn->l = sso->l;
00336        ssn->r = sso->l = hs;
00337        return(ssn);
00338 }
00339 
00340 
00341 static struct source *
00342 splitsource(                /* divide source in two if it's big and long */
00343        struct source *so
00344 )
00345 {
00346        LRSUM  lr;
00347        LRLIN  fit;
00348        register struct srcspan     *ss, *ssn;
00349        struct srcspan       *ssl, *ssnl, head;
00350        int    h;
00351        double mh, mv;
00352        struct source *sn;
00353 
00354        lrclear(&lr);
00355        for (ss = so->first; ss != NULL; ss = ss->next)
00356               for (h = ss->l; h < ss->r; h++)
00357                      lrpoint(h, ss->v, &lr);
00358        if ((double)lr.n/(sampdens*sampdens) < SABIG)
00359               return(NULL);               /* too small */
00360        if (lrfit(&fit, &lr) < 0)
00361               return(NULL);               /* can't fit a line */
00362        if (fit.correlation < LCORR && fit.correlation > -LCORR)
00363               return(NULL);
00364        if (verbose)
00365               fprintf(stderr, "%s: splitting large source\n", progname);
00366        mh = lrxavg(&lr);
00367        mv = lryavg(&lr);
00368        sn = (struct source *)malloc(sizeof(struct source));
00369        if (sn == NULL)
00370               memerr("source records in splitsource");
00371        sn->dom = 0.0;
00372        sn->first = NULL;
00373        ssnl = NULL;
00374        head.next = so->first;
00375        ssl = &head;
00376        for (ss = so->first; ss != NULL; ssl = ss, ss = ss->next)
00377               if ((ssn = splitspan(ss, mh, mv, fit.slope)) != NULL) {
00378                      if (ssn == ss) {     /* remove from old */
00379                             ssl->next = ss->next;
00380                             ss = ssl;
00381                      }
00382                      if (ssnl == NULL)    /* add to new */
00383                             sn->first = ssn;
00384                      else
00385                             ssnl->next = ssn;
00386                      ssn->next = NULL;
00387                      ssnl = ssn;
00388               }
00389        so->first = head.next;
00390        return(sn);
00391 }
00392 
00393 
00394 static void
00395 donesource(                 /* finished with this source */
00396        register struct source      *sp
00397 )
00398 {
00399        struct source *newsrc;
00400        register struct srcspan     *ss;
00401        int    h, n;
00402        double hsum, vsum, d;
00403 
00404        while ((newsrc = splitsource(sp)) != NULL)       /* split it? */
00405               donesource(newsrc);
00406        sp->dom = 0.0;
00407        hsum = vsum = 0.0;
00408        sp->brt = 0.0;
00409        n = 0;
00410        for (ss = sp->first; ss != NULL; ss = ss->next) {
00411               sp->brt += ss->brsum;
00412               n += ss->r - ss->l;
00413               for (h = ss->l; h < ss->r; h++) {
00414                      d = pixsize(h, ss->v);
00415                      hsum += d*h;
00416                      vsum += d*ss->v;
00417                      sp->dom += d;
00418               }
00419        }
00420        freespans(sp);
00421        if (sp->dom <= FTINY) {            /* must be right at edge of image */
00422               free((void *)sp);
00423               return;
00424        }
00425        sp->brt /= (double)n;
00426        compdir(sp->dir, (int)(hsum/sp->dom), (int)(vsum/sp->dom));
00427        sp->next = donelist;
00428        donelist = sp;
00429        if (verbose)
00430               fprintf(stderr,
00431        "%s: source at [%.3f,%.3f,%.3f], dw %.5f, br %.1f (%d samps)\n",
00432                      progname, sp->dir[0], sp->dir[1], sp->dir[2], 
00433                      sp->dom, sp->brt, n);
00434 }
00435 
00436 
00437 static struct source *
00438 findbuddy(                  /* find close enough source to s in l*/
00439        register struct source      *s,
00440        register struct source      *l
00441 )
00442 {
00443        struct source *bestbuddy = NULL;
00444        double d, r, mindist = MAXBUDDY;
00445 
00446        r = sqrt(s->dom/PI);
00447        for ( ; l != NULL; l = l->next) {
00448               d = sqrt(dist2(l->dir, s->dir)) - sqrt(l->dom/PI) - r;
00449               if (d < mindist) {
00450                      bestbuddy = l;
00451                      mindist = d;
00452               }
00453        }
00454        return(bestbuddy);
00455 }
00456 
00457 
00458 extern void
00459 absorb_specks(void)                /* eliminate too-small sources */
00460 {
00461        struct source head, *buddy;
00462        register struct source      *last, *this;
00463 
00464        if (verbose)
00465               fprintf(stderr, "%s: absorbing small sources...\n", progname);
00466        head.next = donelist;
00467        last = &head;
00468        for (this = head.next; this != NULL; this = this->next)
00469               if (TOOSMALL(this)) {
00470                      last->next = this->next;
00471                      buddy = findbuddy(this, head.next);
00472                      if (buddy != NULL)
00473                             mergesource(buddy, this);
00474                      else
00475                             absorb(this);
00476                      this = last;
00477               } else
00478                      last = this;
00479        donelist = head.next;
00480 }
00481 
00482 
00483 static void
00484 absorb(                     /* absorb a source into indirect */
00485        register struct source      *s
00486 )
00487 {
00488        FVECT  dir;
00489        double d;
00490        register int  i;
00491 
00492        for (i = 0; i < nglardirs; i++) {
00493               spinvector(dir, ourview.vdir, ourview.vup, indirect[i].theta);
00494               d = DOT(dir,s->dir)*s->dom*(sampdens*sampdens);
00495               if (d <= 0.0)
00496                      continue;
00497               indirect[i].sum += d * s->brt;
00498               indirect[i].n += d;
00499        }
00500        freespans(s);
00501        free((void *)s);
00502 }
00503 
00504 
00505 static void
00506 freespans(                  /* free spans associated with source */
00507        struct source *sp
00508 )
00509 {
00510        register struct srcspan     *ss;
00511 
00512        while ((ss = sp->first) != NULL) {
00513               sp->first = ss->next;
00514               free((void *)ss);
00515        }
00516 }