--- ray/src/hd/rhd_qtree.c 1997/11/20 11:38:26 3.2 +++ ray/src/hd/rhd_qtree.c 1997/12/05 09:40:05 3.10 @@ -10,25 +10,32 @@ static char SCCSid[] = "$SunId$ SGI"; #include "standard.h" #include "rhd_qtree.h" + /* quantity of leaves to free at a time */ +#ifndef LFREEPCT +#define LFREEPCT 25 +#endif + /* maximum allowed angle difference (deg.) */ +#ifndef MAXANG +#define MAXANG 20. +#endif +#define MAXDIFF2 (long)( MAXANG*MAXANG /90./90.*(1L<<15)*(1L<<15)) + +#define abs(i) ((i) < 0 ? -(i) : (i)) + RTREE qtrunk; /* our quadtree trunk */ -double qtDepthEps = .02; /* epsilon to compare depths (z fraction) */ +double qtDepthEps = .05; /* epsilon to compare depths (z fraction) */ int qtMinNodesiz = 2; /* minimum node dimension (pixels) */ +struct rleaves qtL; /* our pile of leaves */ -static RLEAF *leafpile; /* our collection of leaf values */ -static int nleaves; /* count of leaves in our pile */ -static int bleaf, tleaf; /* bottom and top (next) leaf index (ring) */ - #define TBUNDLESIZ 409 /* number of twigs in a bundle */ static RTREE **twigbundle; /* free twig blocks (NULL term.) */ static int nexttwig; /* next free twig */ -static RTREE emptytree; /* empty tree for test below */ +#define ungetleaf(li) (qtL.tl=(li)) /* dangerous if used improperly */ -#define is_stump(t) (!bcmp((char *)(t), (char *)&emptytree, sizeof(RTREE))) - static RTREE * newtwig() /* allocate a twig */ { @@ -58,20 +65,18 @@ memerr: } -static -freetree(really) /* free allocated twigs */ +qtFreeTree(really) /* free allocated twigs */ int really; { register int i; - if (tmTop != NULL) - tmClearHisto(); - bzero((char *)&qtrunk, sizeof(RTREE)); - nexttwig = 0; + qtrunk.flgs = CH_ANY; /* chop down tree */ if (twigbundle == NULL) return; + i = (TBUNDLESIZ-1+nexttwig)/TBUNDLESIZ; + nexttwig = 0; if (!really) { /* just clear allocated blocks */ - for (i = 0; twigbundle[i] != NULL; i++) + while (i--) bzero((char *)twigbundle[i], TBUNDLESIZ*sizeof(RTREE)); return; } @@ -83,53 +88,67 @@ int really; } -static RLEAF * +static int newleaf() /* allocate a leaf from our pile */ { - if (tleaf++ >= nleaves) /* get next leaf in ring */ - tleaf = 0; - if (tleaf == bleaf) /* need to shake some free */ + int li; + + li = qtL.tl++; + if (qtL.tl >= qtL.nl) /* get next leaf in ring */ + qtL.tl = 0; + if (qtL.tl == qtL.bl) /* need to shake some free */ qtCompost(LFREEPCT); - return(leafpile + tleaf); + return(li); } +#define LEAFSIZ (3*sizeof(float)+2*sizeof(short)+\ + sizeof(TMbright)+6*sizeof(BYTE)) + int qtAllocLeaves(n) /* allocate space for n leaves */ -int n; +register int n; { unsigned nbytes; register unsigned i; - freetree(0); /* make sure tree is empty */ + qtFreeTree(0); /* make sure tree is empty */ if (n <= 0) return(0); - if (nleaves >= n) - return(nleaves); - else if (nleaves > 0) - free((char *)leafpile); + if (qtL.nl >= n) + return(qtL.nl); + else if (qtL.nl > 0) + free(qtL.base); /* round space up to nearest power of 2 */ - nbytes = n*sizeof(RLEAF) + 8; + nbytes = n*LEAFSIZ + 8; for (i = 1024; nbytes > i; i <<= 1) ; - n = (i - 8) / sizeof(RLEAF); - leafpile = (RLEAF *)malloc(n*sizeof(RLEAF)); - if (leafpile == NULL) - return(-1); - nleaves = n; - bleaf = tleaf = 0; - return(nleaves); + n = (i - 8) / LEAFSIZ; /* should we make sure n is even? */ + qtL.base = (char *)malloc(n*LEAFSIZ); + if (qtL.base == NULL) + return(0); + /* assign larger alignment types earlier */ + qtL.wp = (float (*)[3])qtL.base; + qtL.wd = (short (*)[2])(qtL.wp + n); + qtL.brt = (TMbright *)(qtL.wd + n); + qtL.chr = (BYTE (*)[3])(qtL.brt + n); + qtL.rgb = (BYTE (*)[3])(qtL.chr + n); + qtL.nl = n; + qtL.tml = qtL.bl = qtL.tl = 0; + return(n); } +#undef LEAFSIZ + qtFreeLeaves() /* free our allocated leaves and twigs */ { - freetree(1); /* free tree also */ - if (nleaves <= 0) + qtFreeTree(1); /* free tree also */ + if (qtL.nl <= 0) return; - free((char *)leafpile); - leafpile = NULL; - nleaves = 0; + free(qtL.base); + qtL.base = NULL; + qtL.nl = 0; } @@ -140,15 +159,16 @@ register RTREE *tp; register int i, li; for (i = 0; i < 4; i++) - if (tp->flgs & BRF(i)) + if (tp->flgs & BRF(i)) { shaketree(tp->k[i].b); - else if (tp->k[i].l != NULL) { - li = tp->k[i].l - leafpile; - if (bleaf < tleaf ? (li < bleaf || li >= tleaf) : - (li < bleaf && li >= tleaf)) { - tmAddHisto(&tp->k[i].l->brt, 1, -1); - tp->k[i].l = NULL; - } + if (is_stump(tp->k[i].b)) + tp->flgs &= ~BRF(i); + } else if (tp->flgs & LFF(i)) { + li = tp->k[i].li; + if (qtL.bl < qtL.tl ? + (li < qtL.bl || li >= qtL.tl) : + (li < qtL.bl && li >= qtL.tl)) + tp->flgs &= ~LFF(i); } } @@ -157,41 +177,157 @@ int qtCompost(pct) /* free up some leaves */ int pct; { - int nused, nclear; + int nused, nclear, nmapped; + /* figure out how many leaves to clear */ - nclear = nleaves * pct / 100; + nclear = qtL.nl * pct / 100; + nused = qtL.tl - qtL.bl; + if (nused <= 0) nused += qtL.nl; + nclear -= qtL.nl - nused; if (nclear <= 0) return(0); - nused = tleaf > bleaf ? tleaf-bleaf : tleaf+nleaves-bleaf; if (nclear >= nused) { /* clear them all */ - freetree(0); - bleaf = tleaf = 0; + qtFreeTree(0); + qtL.tml = qtL.bl = qtL.tl = 0; return(nused); } /* else clear leaves from bottom */ - bleaf = (bleaf + nclear) % nleaves; + nmapped = qtL.tml - qtL.bl; + if (nmapped < 0) nmapped += qtL.nl; + qtL.bl += nclear; + if (qtL.bl >= qtL.nl) qtL.bl -= qtL.nl; + if (nmapped <= nclear) qtL.tml = qtL.bl; shaketree(&qtrunk); return(nclear); } static -addleaf(lp) /* add a leaf to our tree */ -RLEAF *lp; +encodedir(pa, dv) /* encode a normalized direction vector */ +short pa[2]; +FVECT dv; { + pa[1] = 0; + if (dv[2] >= 1.) + pa[0] = (1L<<15)-1; + else if (dv[2] <= -1.) + pa[0] = -((1L<<15)-1); + else { + pa[0] = ((1L<<15)-1)/(PI/2.) * asin(dv[2]); + pa[1] = ((1L<<15)-1)/PI * atan2(dv[1], dv[0]); + } +} + + +#define ALTSHFT 5 +#define NALT (1<>(15-ALTSHFT)] + +static unsigned short azisftab[NALT]; + +static +azisftinit(alt) /* initialize azimuth scale factor table */ +int alt; +{ + register int i; + + for (i = NALT; i--; ) + azisftab[i] = 2.*(1L<<15) * cos(PI/2.*(i+.5)/NALT); + return(azisft(alt)); +} + +#define azisf(alt) (azisftab[0] ? azisft(alt) : azisftinit(alt)) >> 15 + +static long +dir2diff(pa1, pa2) /* relative distance^2 between directions */ +short pa1[2], pa2[2]; +{ + long altd2, azid2; + int alt; + + altd2 = pa1[0] - pa2[0]; /* get altitude difference^2 */ + altd2 *= altd2; + if (altd2 > MAXDIFF2) + return(altd2); /* too large already */ + azid2 = pa1[1] - pa2[1]; /* get adjusted azimuth difference^2 */ + if (azid2 < 0) azid2 = -azid2; + if (azid2 >= 1L<<15) { /* wrap sphere */ + azid2 -= 1L<<16; + if (azid2 < 0) azid2 = -azid2; + } + alt = (pa1[0]+pa2[0])/2; + azid2 = azid2*azisf(alt); /* evaluation order is important */ + azid2 *= azid2; + return(altd2 + azid2); +} + + +int +qtFindLeaf(x, y) /* find closest leaf to (x,y) */ +int x, y; +{ register RTREE *tp = &qtrunk; + int li = -1; int x0=0, y0=0, x1=odev.hres, y1=odev.vres; - RLEAF *lo = NULL; + int mx, my; + register int q; + /* check limits */ + if (x < 0 || x >= odev.hres || y < 0 || y >= odev.vres) + return(-1); + /* find nearby leaf in our tree */ + for ( ; ; ) { + for (q = 0; q < 4; q++) /* find any leaf this level */ + if (tp->flgs & LFF(q)) { + li = tp->k[q].li; + break; + } + q = 0; /* which quadrant are we? */ + mx = (x0 + x1) >> 1; + my = (y0 + y1) >> 1; + if (x < mx) x1 = mx; + else {x0 = mx; q |= 01;} + if (y < my) y1 = my; + else {y0 = my; q |= 02;} + if (tp->flgs & BRF(q)) { /* branch down if not a leaf */ + tp = tp->k[q].b; + continue; + } + if (tp->flgs & LFF(q)) /* good shot! */ + return(tp->k[q].li); + return(li); /* else return what we have */ + } +} + + +static +addleaf(li) /* add a leaf to our tree */ +int li; +{ + register RTREE *tp = &qtrunk; + int x0=0, y0=0, x1=odev.hres, y1=odev.vres; + int lo = -1; + long d2; + short dc[2]; int x, y, mx, my; double z; FVECT ip, wp; register int q; - /* compute leaf location */ - VCOPY(wp, lp->wp); + /* compute leaf location in view */ + VCOPY(wp, qtL.wp[li]); viewloc(ip, &odev.v, wp); if (ip[2] <= 0. || ip[0] < 0. || ip[0] >= 1. || ip[1] < 0. || ip[1] >= 1.) - return; + return(0); /* behind or outside view */ +#ifdef DEBUG + if (odev.v.type == VT_PAR | odev.v.vfore > FTINY) + error(INTERNAL, "bad view assumption in addleaf"); +#endif + for (q = 0; q < 3; q++) + wp[q] = (wp[q] - odev.v.vp[q])/ip[2]; + encodedir(dc, wp); /* compute pixel direction */ + d2 = dir2diff(dc, qtL.wd[li]); + if (d2 > MAXDIFF2) + return(0); /* leaf dir. too far off */ x = ip[0] * odev.hres; y = ip[1] * odev.vres; z = ip[2]; @@ -209,203 +345,107 @@ RLEAF *lp; tp = tp->k[q].b; continue; } - if (tp->k[q].l == NULL) { /* found stem for leaf */ - tp->k[q].l = lp; - tp->flgs |= CHF(q); + if (!(tp->flgs & LFF(q))) { /* found stem for leaf */ + tp->k[q].li = li; + tp->flgs |= CHLFF(q); break; } - /* check existing leaf */ - if (lo != tp->k[q].l) { - lo = tp->k[q].l; - VCOPY(wp, lo->wp); + if (lo != tp->k[q].li) { /* check old leaf */ + lo = tp->k[q].li; + VCOPY(wp, qtL.wp[lo]); viewloc(ip, &odev.v, wp); } /* is node minimum size? */ - if (x1-x0 <= qtMinNodesiz || y1-y0 <= qtMinNodesiz) { - if (z > (1.-qtDepthEps)*ip[2]) /* who is closer? */ - return; /* old one is */ - tp->k[q].l = lp; /* new one is */ + if (y1-y0 <= qtMinNodesiz || x1-x0 <= qtMinNodesiz) { + if (z > (1.+qtDepthEps)*ip[2]) + return(0); /* old one closer */ + if (z >= (1.-qtDepthEps)*ip[2] && + dir2diff(dc, qtL.wd[lo]) < d2) + return(0); /* old one better */ + tp->k[q].li = li; /* else new one is */ tp->flgs |= CHF(q); - tmAddHisto(&lo->brt, 1, -1); /* drop old one */ break; } - tp->flgs |= CHBRF(q); /* else grow tree */ + tp->flgs &= ~LFF(q); /* else grow tree */ + tp->flgs |= CHBRF(q); tp = tp->k[q].b = newtwig(); - tp->flgs |= CH_ANY; /* all new */ q = 0; /* old leaf -> new branch */ mx = ip[0] * odev.hres; my = ip[1] * odev.vres; if (mx >= (x0 + x1) >> 1) q |= 01; if (my >= (y0 + y1) >> 1) q |= 02; - tp->k[q].l = lo; + tp->k[q].li = lo; + tp->flgs |= LFF(q)|CH_ANY; /* all new */ } - tmAddHisto(&lp->brt, 1, 1); /* add leaf to histogram */ + return(1); /* done */ } -dev_value(c, p) /* add a pixel value to our output queue */ +dev_value(c, p, v) /* add a pixel value to our quadtree */ COLR c; -FVECT p; +FVECT p, v; { - register RLEAF *lp; + register int li; - lp = newleaf(); - VCOPY(lp->wp, p); - tmCvColrs(&lp->brt, lp->chr, c, 1); - addleaf(lp); + li = newleaf(); + VCOPY(qtL.wp[li], p); + encodedir(qtL.wd[li], v); + tmCvColrs(&qtL.brt[li], qtL.chr[li], c, 1); + if (!addleaf(li)) + ungetleaf(li); } qtReplant() /* replant our tree using new view */ { register int i; - - if (bleaf == tleaf) /* anything to replant? */ + /* anything to replant? */ + if (qtL.bl == qtL.tl) return; - freetree(0); /* blow the tree away */ - /* now rebuild it */ - for (i = bleaf; i != tleaf; ) { - addleaf(leafpile+i); - if (++i >= nleaves) i = 0; + qtFreeTree(0); /* blow the old tree away */ + /* regrow it in new place */ + for (i = qtL.bl; i != qtL.tl; ) { + addleaf(i); + if (++i >= qtL.nl) i = 0; } - tmComputeMapping(0., 0., 0.); /* update the display */ - qtUpdate(); } -static -redraw(ca, tp, x0, y0, x1, y1, l) /* redraw portion of a tree */ -BYTE ca[3]; /* returned average color */ -register RTREE *tp; -int x0, y0, x1, y1; -int l[2][2]; +qtMapLeaves(redo) /* map our leaves to RGB */ +int redo; { - int csm[3], nc; - BYTE rgb[3]; - int quads = CH_ANY; - int mx, my; - register int i; - /* compute midpoint */ - mx = (x0 + x1) >> 1; - my = (y0 + y1) >> 1; - /* see what to do */ - if (l[0][0] >= mx) - quads &= ~(CHF(2)|CHF(0)); - else if (l[0][1] <= mx) - quads &= ~(CHF(3)|CHF(1)); - if (l[1][0] >= my) - quads &= ~(CHF(1)|CHF(0)); - else if (l[1][1] <= my) - quads &= ~(CHF(3)|CHF(2)); - tp->flgs &= ~quads; /* mark them done */ - csm[0] = csm[1] = csm[2] = nc = 0; - /* do leaves first */ - for (i = 0; i < 4; i++) - if (quads & CHF(i) && !(tp->flgs & BRF(i)) && - tp->k[i].l != NULL) { - tmMapPixels(rgb, &tp->k[i].l->brt, tp->k[i].l->chr, 1); - dev_paintr(rgb, i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my); - csm[0] += rgb[0]; csm[1] += rgb[1]; csm[2] += rgb[2]; - nc++; - quads &= ~CHF(i); - } - /* now do branches */ - for (i = 0; i < 4; i++) - if (quads & CHF(i) && tp->flgs & BRF(i)) { - redraw(rgb, tp->k[i].b, i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my, l); - csm[0] += rgb[0]; csm[1] += rgb[1]; csm[2] += rgb[2]; - nc++; - quads &= ~CHF(i); - } - if (nc > 1) { - ca[0] = csm[0]/nc; ca[1] = csm[1]/nc; ca[2] = csm[2]/nc; + int aorg, alen, borg, blen; + /* recompute mapping? */ + if (redo) + qtL.tml = qtL.bl; + /* already done? */ + if (qtL.tml == qtL.tl) + return(1); + /* compute segments */ + aorg = qtL.tml; + if (qtL.tl >= aorg) { + alen = qtL.tl - aorg; + blen = 0; } else { - ca[0] = csm[0]; ca[1] = csm[1]; ca[2] = csm[2]; + alen = qtL.nl - aorg; + borg = 0; + blen = qtL.tl; } - if (!quads) return; - /* fill in gaps with average */ - for (i = 0; i < 4; i++) - if (quads & CHF(i)) - dev_paintr(ca, i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my); -} - - -static -update(ca, tp, x0, y0, x1, y1) /* update tree display as needed */ -BYTE ca[3]; /* returned average color */ -register RTREE *tp; -int x0, y0, x1, y1; -{ - int csm[3], nc; - BYTE rgb[3]; - int gaps = 0; - int mx, my; - register int i; - /* compute midpoint */ - mx = (x0 + x1) >> 1; - my = (y0 + y1) >> 1; - csm[0] = csm[1] = csm[2] = nc = 0; - /* do leaves first */ - for (i = 0; i < 4; i++) - if ((tp->flgs & CHBRF(i)) == CHF(i)) { - if (tp->k[i].l == NULL) { - gaps |= 1<k[i].l->brt, tp->k[i].l->chr, 1); - dev_paintr(rgb, i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my); - csm[0] += rgb[0]; csm[1] += rgb[1]; csm[2] += rgb[2]; - nc++; - } - /* now do branches */ - for (i = 0; i < 4; i++) - if ((tp->flgs & CHBRF(i)) == CHBRF(i)) { - update(rgb, tp->k[i].b, i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my); - csm[0] += rgb[0]; csm[1] += rgb[1]; csm[2] += rgb[2]; - nc++; - } - if (nc > 1) { - ca[0] = csm[0]/nc; ca[1] = csm[1]/nc; ca[2] = csm[2]/nc; - } else { - ca[0] = csm[0]; ca[1] = csm[1]; ca[2] = csm[2]; + /* (re)compute tone mapping? */ + if (qtL.tml == qtL.bl) { + tmClearHisto(); + tmAddHisto(qtL.brt+aorg, alen, 1); + if (blen > 0) + tmAddHisto(qtL.brt+borg, blen, 1); + if (tmComputeMapping(0., 0., 0.) != TM_E_OK) + return(0); } - /* fill in gaps with average */ - for (i = 0; gaps && i < 4; gaps >>= 1, i++) - if (gaps & 01) - dev_paintr(ca, i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my); - tp->flgs &= ~CH_ANY; /* all done */ -} - - -qtRedraw(x0, y0, x1, y1) /* redraw part of our screen */ -int x0, y0, x1, y1; -{ - int lim[2][2]; - BYTE ca[3]; - - if (is_stump(&qtrunk)) - return; - if ((lim[0][0]=x0) == 0 & (lim[1][0]=y0) == 0 & - (lim[0][1]=x1) == odev.hres & (lim[1][1]=y1) == odev.vres || - tmTop->lumap == NULL) - tmComputeMapping(0., 0., 0.); - redraw(ca, &qtrunk, 0, 0, odev.hres, odev.vres, lim); -} - - -qtUpdate() /* update our tree display */ -{ - BYTE ca[3]; - - if (is_stump(&qtrunk)) - return; - if (tmTop->lumap == NULL) - tmComputeMapping(0., 0., 0.); - update(ca, &qtrunk, 0, 0, odev.hres, odev.vres); + if (tmMapPixels(qtL.rgb+aorg, qtL.brt+aorg, + qtL.chr+aorg, alen) != TM_E_OK) + return(0); + if (blen > 0) + tmMapPixels(qtL.rgb+borg, qtL.brt+borg, + qtL.chr+borg, blen); + qtL.tml = qtL.tl; + return(1); }