--- ray/src/hd/rhd_qtree.c	1997/11/25 16:52:04	3.9
+++ 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,9 @@ struct rleaves	qtL;		/* our pile of leaves */
 static RTREE	**twigbundle;	/* free twig blocks (NULL term.) */
 static int	nexttwig;	/* next free twig */
 
+#define ungetleaf(li)	(qtL.tl=(li))	/* dangerous if used improperly */
 
+
 static RTREE *
 newtwig()			/* allocate a twig */
 {
@@ -92,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 */
@@ -118,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;
@@ -190,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<<ALTSHFT)
+#define azisft(alt)	azisftab[abs(alt)>>(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;
@@ -234,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];
@@ -266,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;
 		}
@@ -291,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 quadtree */
+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);
 }