--- ray/src/hd/rhd_qtree.c 1997/11/25 11:15:20 3.7 +++ ray/src/hd/rhd_qtree.c 1997/12/05 09:40:05 3.10 @@ -14,9 +14,17 @@ static char SCCSid[] = "$SunId$ SGI"; #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 */ @@ -25,7 +33,7 @@ struct rleaves qtL; /* our pile of leaves */ static RTREE **twigbundle; /* free twig blocks (NULL term.) */ static int nexttwig; /* next free twig */ -#define is_stump(t) (!((t)->flgs & (BR_ANY|LF_ANY))) +#define ungetleaf(li) (qtL.tl=(li)) /* dangerous if used improperly */ static RTREE * @@ -94,7 +102,8 @@ newleaf() /* allocate a leaf from our pile */ } -#define LEAFSIZ (3*sizeof(float)+sizeof(TMbright)+6*sizeof(BYTE)) +#define LEAFSIZ (3*sizeof(float)+2*sizeof(short)+\ + sizeof(TMbright)+6*sizeof(BYTE)) int qtAllocLeaves(n) /* allocate space for n leaves */ @@ -120,7 +129,8 @@ register int n; return(0); /* assign larger alignment types earlier */ qtL.wp = (float (*)[3])qtL.base; - qtL.brt = (TMbright *)(qtL.wp + n); + 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; @@ -192,6 +202,66 @@ int pct; } +static +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; @@ -236,16 +306,28 @@ 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 */ + /* 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]; @@ -268,17 +350,19 @@ int li; tp->flgs |= CHLFF(q); break; } - /* check existing leaf */ - if (lo != tp->k[q].li) { + 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].li = li; /* 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); break; } @@ -293,19 +377,22 @@ int li; tp->k[q].li = lo; tp->flgs |= LFF(q)|CH_ANY; /* all new */ } + 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 int li; li = newleaf(); VCOPY(qtL.wp[li], p); + encodedir(qtL.wd[li], v); tmCvColrs(&qtL.brt[li], qtL.chr[li], c, 1); - addleaf(li); + if (!addleaf(li)) + ungetleaf(li); } @@ -361,111 +448,4 @@ int redo; qtL.chr+borg, blen); qtL.tml = qtL.tl; return(1); -} - - -static -redraw(tp, x0, y0, x1, y1, l) /* mark portion of a tree for redraw */ -register RTREE *tp; -int x0, y0, x1, y1; -int l[2][2]; -{ - 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 quadrants for update */ - /* climb the branches */ - for (i = 0; i < 4; i++) - if (tp->flgs & BRF(i) && quads & CHF(i)) - redraw(tp->k[i].b, i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my, l); -} - - -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; - register BYTE *cp; - 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 & CHF(i))) - continue; - if (tp->flgs & LFF(i)) { - dev_paintr(cp=qtL.rgb[tp->k[i].li], - i&01 ? mx : x0, i&02 ? my : y0, - i&01 ? x1 : mx, i&02 ? y1 : my); - csm[0] += cp[0]; csm[1] += cp[1]; csm[2] += cp[2]; - nc++; - } else if (!(tp->flgs & BRF(i))) - gaps |= 1<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]; - } - /* 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 or all of our screen */ -int x0, y0, x1, y1; -{ - int lim[2][2]; - - if (is_stump(&qtrunk)) - return; - if (!qtMapLeaves((lim[0][0]=x0) <= 0 & (lim[1][0]=y0) <= 0 & - (lim[0][1]=x1) >= odev.hres-1 & (lim[1][1]=y1) >= odev.vres-1)) - return; - redraw(&qtrunk, 0, 0, odev.hres, odev.vres, lim); -} - - -qtUpdate() /* update our tree display */ -{ - BYTE ca[3]; - - if (is_stump(&qtrunk)) - return; - if (!qtMapLeaves(0)) - return; - update(ca, &qtrunk, 0, 0, odev.hres, odev.vres); }