Back to index

radiance  4R0+20100331
rhd_qtree.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: rhd_qtree.c,v 3.26 2005/01/07 20:33:02 greg Exp $";
00003 #endif
00004 /*
00005  * Quadtree driver support routines.
00006  */
00007 
00008 #include <string.h>
00009 
00010 #include "standard.h"
00011 #include "rhd_qtree.h"
00012                             /* quantity of leaves to free at a time */
00013 #ifndef LFREEPCT
00014 #define       LFREEPCT      25
00015 #endif
00016                             /* maximum allowed angle difference (deg.) */
00017 #ifndef MAXANG
00018 #define MAXANG              20
00019 #endif
00020 #if MAXANG>0
00021 #define MAXDIFF2     ( MAXANG*MAXANG * (PI*PI/180./180.))
00022 #endif
00023 
00024 #define abs(i)              ((i) < 0 ? -(i) : (i))
00025 
00026 RTREE  qtrunk;                     /* our quadtree trunk */
00027 double qtDepthEps = .05;    /* epsilon to compare depths (z fraction) */
00028 int    qtMinNodesiz = 2;    /* minimum node dimension (pixels) */
00029 struct rleaves       qtL;          /* our pile of leaves */
00030 
00031 int    rayqleft = 0;        /* rays left to queue before flush */
00032 
00033 static int32  falleaves;    /* our list of fallen leaves */
00034 
00035 #define composted(li)       (qtL.bl <= qtL.tl ? \
00036                                    ((li) < qtL.bl || (li) >= qtL.tl) : \
00037                                    ((li) < qtL.bl && (li) >= qtL.tl))
00038 
00039 #define       TBUNDLESIZ    409    /* number of twigs in a bundle */
00040 
00041 static RTREE  **twigbundle; /* free twig blocks (NULL term.) */
00042 static int    nexttwig;     /* next free twig */
00043 
00044 static RTREE *newtwig(void);
00045 static void qtFreeTree(int really);
00046 static void shaketree(RTREE *tp);
00047 static int putleaf(int li, int drop);
00048 
00049 
00050 static RTREE *
00051 newtwig(void)               /* allocate a twig */
00052 {
00053        register int  bi;
00054 
00055        if (twigbundle == NULL) {   /* initialize */
00056               twigbundle = (RTREE **)malloc(sizeof(RTREE *));
00057               if (twigbundle == NULL)
00058                      goto memerr;
00059               twigbundle[0] = NULL;
00060        }
00061        bi = nexttwig / TBUNDLESIZ;
00062        if (twigbundle[bi] == NULL) {      /* new block */
00063               twigbundle = (RTREE **)realloc((void *)twigbundle,
00064                                    (bi+2)*sizeof(RTREE *));
00065               if (twigbundle == NULL)
00066                      goto memerr;
00067               twigbundle[bi] = (RTREE *)calloc(TBUNDLESIZ, sizeof(RTREE));
00068               if (twigbundle[bi] == NULL)
00069                      goto memerr;
00070               twigbundle[bi+1] = NULL;
00071        }
00072                             /* nexttwig++ % TBUNDLESIZ */
00073        return(twigbundle[bi] + (nexttwig++ - bi*TBUNDLESIZ));
00074 memerr:
00075        error(SYSTEM, "out of memory in newtwig");
00076        return NULL; /* pro forma return */
00077 }
00078 
00079 
00080 static void
00081 qtFreeTree(          /* free allocated twigs */
00082        int    really
00083 )
00084 {
00085        register int  i;
00086 
00087        qtrunk.flgs = CH_ANY;       /* chop down tree */
00088        if (twigbundle == NULL)
00089               return;
00090        i = (TBUNDLESIZ-1+nexttwig)/TBUNDLESIZ;
00091        nexttwig = 0;
00092        if (!really) {              /* just clear allocated blocks */
00093               while (i--)
00094                      memset((char *)twigbundle[i], '\0', TBUNDLESIZ*sizeof(RTREE));
00095               return;
00096        }
00097                             /* else "really" means free up memory */
00098        for (i = 0; twigbundle[i] != NULL; i++)
00099               free((void *)twigbundle[i]);
00100        free((void *)twigbundle);
00101        twigbundle = NULL;
00102 }
00103 
00104 
00105 #define       LEAFSIZ              (3*sizeof(float)+sizeof(int32)+\
00106                      sizeof(TMbright)+6*sizeof(BYTE))
00107 
00108 extern int
00109 qtAllocLeaves(              /* allocate space for n leaves */
00110        register int  n
00111 )
00112 {
00113        unsigned      nbytes;
00114        register unsigned    i;
00115 
00116        qtFreeTree(0);              /* make sure tree is empty */
00117        if (n <= 0)
00118               return(0);
00119        if (qtL.nl >= n)
00120               return(qtL.nl);
00121        else if (qtL.nl > 0)
00122               free(qtL.base);
00123                             /* round space up to nearest power of 2 */
00124        nbytes = n*LEAFSIZ + 8;
00125        for (i = 1024; nbytes > i; i <<= 1)
00126               ;
00127        n = (i - 8) / LEAFSIZ;      /* should we make sure n is even? */
00128        qtL.base = (char *)malloc(n*LEAFSIZ);
00129        if (qtL.base == NULL)
00130               return(0);
00131                             /* assign larger alignment types earlier */
00132        qtL.wp = (float (*)[3])qtL.base;
00133        qtL.wd = (int32 *)(qtL.wp + n);
00134        qtL.brt = (TMbright *)(qtL.wd + n);
00135        qtL.chr = (BYTE (*)[3])(qtL.brt + n);
00136        qtL.rgb = (BYTE (*)[3])(qtL.chr + n);
00137        qtL.nl = n;
00138        qtL.tml = qtL.bl = qtL.tl = 0;
00139        falleaves = -1;
00140        return(n);
00141 }
00142 
00143 #undef LEAFSIZ
00144 
00145 
00146 extern void
00147 qtFreeLeaves(void)                 /* free our allocated leaves and twigs */
00148 {
00149        qtFreeTree(1);              /* free tree also */
00150        if (qtL.nl <= 0)
00151               return;
00152        free(qtL.base);
00153        qtL.base = NULL;
00154        qtL.nl = 0;
00155 }
00156 
00157 
00158 static void
00159 shaketree(                  /* shake dead leaves from tree */
00160        register RTREE       *tp
00161 )
00162 {
00163        register int  i, li;
00164 
00165        for (i = 0; i < 4; i++)
00166               if (tp->flgs & BRF(i)) {
00167                      shaketree(tp->k[i].b);
00168                      if (is_stump(tp->k[i].b))
00169                             tp->flgs &= ~BRF(i);
00170               } else if (tp->flgs & LFF(i)) {
00171                      li = tp->k[i].li;
00172                      if (composted(li))
00173                             tp->flgs &= ~LFF(i);
00174               }
00175 }
00176 
00177 
00178 extern int
00179 qtCompost(                  /* free up some leaves */
00180        int    pct
00181 )
00182 {
00183        register int32       *fl;
00184        int    nused, nclear, nmapped;
00185                             /* figure out how many leaves to clear */
00186        nclear = qtL.nl * pct / 100;
00187        nused = qtL.tl - qtL.bl;
00188        if (nused <= 0) nused += qtL.nl;
00189        nclear -= qtL.nl - nused;
00190        if (nclear <= 0)
00191               return(0);
00192        if (nclear >= nused) {      /* clear them all */
00193               qtFreeTree(0);
00194               qtL.tml = qtL.bl = qtL.tl = 0;
00195               falleaves = -1;
00196               return(nused);
00197        }
00198                             /* else clear leaves from bottom */
00199        nmapped = qtL.tml - qtL.bl;
00200        if (nmapped < 0) nmapped += qtL.nl;
00201        qtL.bl += nclear;
00202        if (qtL.bl >= qtL.nl) qtL.bl -= qtL.nl;
00203        if (nmapped <= nclear) qtL.tml = qtL.bl;
00204        shaketree(&qtrunk);  /* dereference composted leaves */
00205        for (fl = &falleaves; *fl >= 0; fl = qtL.wd + *fl)
00206               while (composted(*fl))
00207                      if ((*fl = qtL.wd[*fl]) < 0)
00208                             return(nclear);
00209        return(nclear);
00210 }
00211 
00212 
00213 extern int
00214 qtFindLeaf(          /* find closest leaf to (x,y) */
00215        int    x,
00216        int    y
00217 )
00218 {
00219        register RTREE       *tp = &qtrunk;
00220        int    li = -1;
00221        int    x0=0, y0=0, x1=odev.hres, y1=odev.vres;
00222        int    mx, my;
00223        register int  q;
00224                                    /* check limits */
00225        if (x < 0 || x >= odev.hres || y < 0 || y >= odev.vres)
00226               return(-1);
00227                                    /* find nearby leaf in our tree */
00228        for ( ; ; ) {
00229               for (q = 0; q < 4; q++)            /* find any leaf this level */
00230                      if (tp->flgs & LFF(q)) {
00231                             li = tp->k[q].li;
00232                             break;
00233                      }
00234               q = 0;                      /* which quadrant are we? */
00235               mx = (x0 + x1) >> 1;
00236               my = (y0 + y1) >> 1;
00237               if (x < mx) x1 = mx;
00238               else {x0 = mx; q |= 01;}
00239               if (y < my) y1 = my;
00240               else {y0 = my; q |= 02;}
00241               if (tp->flgs & BRF(q)) {    /* branch down if not a leaf */
00242                      tp = tp->k[q].b;
00243                      continue;
00244               }
00245               if (tp->flgs & LFF(q))             /* good shot! */
00246                      return(tp->k[q].li);
00247               return(li);                 /* else return what we have */
00248        }
00249 }
00250 
00251 
00252 static int
00253 putleaf(             /* put a leaf in our tree */
00254        register int  li,
00255        int    drop
00256 )
00257 {
00258        register RTREE       *tp = &qtrunk;
00259        int    x0=0, y0=0, x1=odev.hres, y1=odev.vres;
00260        register int  lo = -1;
00261        double d2;
00262        int    x, y, mx, my;
00263        double z;
00264        FVECT  ip, wp, vd;
00265        register int  q;
00266                                    /* check for dead leaf */
00267        if (!qtL.chr[li][1] && !(qtL.chr[li][0] | qtL.chr[li][2]))
00268               return(0);
00269                                    /* compute leaf location in view */
00270        VCOPY(wp, qtL.wp[li]);
00271        viewloc(ip, &odev.v, wp);
00272        if (ip[2] <= 0. || ip[0] < 0. || ip[0] >= 1.
00273                      || ip[1] < 0. || ip[1] >= 1.)
00274               goto dropit;                /* behind or outside view */
00275 #ifdef DEBUG
00276        if (odev.v.type == VT_PAR | odev.v.vfore > FTINY)
00277               error(INTERNAL, "bad view assumption in putleaf");
00278 #endif
00279        for (q = 0; q < 3; q++)
00280               vd[q] = (wp[q] - odev.v.vp[q])/ip[2];
00281        d2 = fdir2diff(qtL.wd[li], vd);
00282 #ifdef MAXDIFF2
00283        if (d2 > MAXDIFF2)
00284               goto dropit;                /* leaf dir. too far off */
00285 #endif
00286        x = ip[0] * odev.hres;
00287        y = ip[1] * odev.vres;
00288        z = ip[2];
00289                                    /* find the place for it */
00290        for ( ; ; ) {
00291               q = 0;                      /* which quadrant? */
00292               mx = (x0 + x1) >> 1;
00293               my = (y0 + y1) >> 1;
00294               if (x < mx) x1 = mx;
00295               else {x0 = mx; q |= 01;}
00296               if (y < my) y1 = my;
00297               else {y0 = my; q |= 02;}
00298               if (tp->flgs & BRF(q)) {    /* move to next branch */
00299                      tp->flgs |= CHF(q);         /* not sure; guess */
00300                      tp = tp->k[q].b;
00301                      continue;
00302               }
00303               if (!(tp->flgs & LFF(q))) { /* found stem for leaf */
00304                      tp->k[q].li = li;
00305                      tp->flgs |= CHLFF(q);
00306                      return(1);
00307               }      
00308               if (lo != tp->k[q].li) {    /* check old leaf */
00309                      lo = tp->k[q].li;
00310                      VCOPY(wp, qtL.wp[lo]);
00311                      viewloc(ip, &odev.v, wp);
00312               }
00313                                           /* is node minimum size? */
00314               if (y1-y0 <= qtMinNodesiz || x1-x0 <= qtMinNodesiz) {
00315                      if (z > (1.+qtDepthEps)*ip[2])
00316                             break;               /* old one closer */
00317                      if (z >= (1.-qtDepthEps)*ip[2] &&
00318                                    fdir2diff(qtL.wd[lo], vd) < d2)
00319                             break;               /* old one better */
00320                      tp->k[q].li = li;           /* attach new */
00321                      tp->flgs |= CHF(q);
00322                      li = lo;                    /* drop old... */
00323                      break;
00324               }
00325               tp->flgs &= ~LFF(q);        /* else grow tree */
00326               tp->flgs |= CHBRF(q);
00327               tp = tp->k[q].b = newtwig();
00328               q = 0;                      /* old leaf -> new branch */
00329               mx = ip[0] * odev.hres;
00330               my = ip[1] * odev.vres;
00331               if (mx >= (x0 + x1) >> 1) q |= 01;
00332               if (my >= (y0 + y1) >> 1) q |= 02;
00333               tp->flgs = CH_ANY|LFF(q);   /* all new */
00334               tp->k[q].li = lo;
00335        }
00336 dropit:
00337        if (drop) {
00338               if (li+1 == (qtL.tl ? qtL.tl : qtL.nl))
00339                      qtL.tl = li;         /* special case */
00340               else {
00341                      qtL.chr[li][0] = qtL.chr[li][1] = qtL.chr[li][2] = 0;
00342                      qtL.wd[li] = falleaves;
00343                      falleaves = li;
00344               }
00345        }
00346        return(li == lo);
00347 }
00348 
00349 
00350 extern void
00351 dev_value(           /* add a pixel value to our quadtree */
00352        COLR   c,
00353        FVECT  d,
00354        FVECT  p
00355 )
00356 {
00357        register int  li;
00358        int    mapit;
00359                             /* grab a leaf */
00360        if (!imm_mode && falleaves >= 0) { /* check for fallen leaves */
00361               li = falleaves;
00362               falleaves = qtL.wd[li];
00363               mapit = qtL.tml <= qtL.tl ?
00364                             (li < qtL.tml || li >= qtL.tl) :
00365                             (li < qtL.tml && li >= qtL.tl) ;
00366        } else {                           /* else allocate new one */
00367               li = qtL.tl++;
00368               if (qtL.tl >= qtL.nl)              /* next leaf in ring */
00369                      qtL.tl = 0;
00370               if (qtL.tl == qtL.bl)              /* need to shake some free */
00371                      qtCompost(LFREEPCT);
00372               mapit = 0;                  /* we'll map it later */
00373        }
00374        if (p == NULL)
00375               VSUM(qtL.wp[li], odev.v.vp, d, FHUGE);
00376        else
00377               VCOPY(qtL.wp[li], p);
00378        qtL.wd[li] = encodedir(d);
00379        tmCvColrs(tmGlobal, &qtL.brt[li], qtL.chr[li], (COLR *)c, 1);
00380        if (putleaf(li, 1)) {
00381               if (mapit)
00382                      tmMapPixels(tmGlobal, (BYTE *)(qtL.rgb+li), qtL.brt+li,
00383                                    (BYTE *)(qtL.chr+li), 1);
00384               if (--rayqleft == 0)
00385                      dev_flush();         /* flush output */
00386        }
00387 }
00388 
00389 
00390 extern void
00391 qtReplant(void)                    /* replant our tree using new view */
00392 {
00393        register int  i;
00394                                    /* anything to replant? */
00395        if (qtL.bl == qtL.tl)
00396               return;
00397        qtFreeTree(0);                     /* blow the old tree away */
00398                                    /* regrow it in new place */
00399        for (i = qtL.bl; i != qtL.tl; ) {
00400               putleaf(i, 0);
00401               if (++i >= qtL.nl) i = 0;
00402        }
00403 }
00404 
00405 
00406 extern int
00407 qtMapLeaves(         /* map our leaves to RGB */
00408        int    redo
00409 )
00410 {
00411        int    aorg, alen, borg, blen;
00412                                    /* recompute mapping? */
00413        if (redo)
00414               qtL.tml = qtL.bl;
00415                                    /* already done? */
00416        if (qtL.tml == qtL.tl)
00417               return(1);
00418                                    /* compute segments */
00419        aorg = qtL.tml;
00420        if (qtL.tl >= aorg) {
00421               alen = qtL.tl - aorg;
00422               blen = 0;
00423        } else {
00424               alen = qtL.nl - aorg;
00425               borg = 0;
00426               blen = qtL.tl;
00427        }
00428                                    /* (re)compute tone mapping? */
00429        if (qtL.tml == qtL.bl) {
00430               tmClearHisto(tmGlobal);
00431               tmAddHisto(tmGlobal, qtL.brt+aorg, alen, 1);
00432               if (blen > 0)
00433                      tmAddHisto(tmGlobal, qtL.brt+borg, blen, 1);
00434               if (tmComputeMapping(tmGlobal, 0., 0., 0.) != TM_E_OK)
00435                      return(0);
00436        }
00437        if (tmMapPixels(tmGlobal, (BYTE *)(qtL.rgb+aorg), qtL.brt+aorg,
00438                      (BYTE *)(qtL.chr+aorg), alen) != TM_E_OK)
00439               return(0);
00440        if (blen > 0)
00441               tmMapPixels(tmGlobal, (BYTE *)(qtL.rgb+borg), qtL.brt+borg,
00442                             (BYTE *)(qtL.chr+borg), blen);
00443        qtL.tml = qtL.tl;
00444        return(1);
00445 }