--- ray/src/hd/rhd_qtree.c	1997/11/25 16:52:04	3.9
+++ ray/src/hd/rhd_qtree.c	1997/12/05 16:22:49	3.12
@@ -14,9 +14,18 @@ static char SCCSid[] = "$SunId$ SGI";
 #ifndef LFREEPCT
 #define	LFREEPCT	25
 #endif
+				/* maximum allowed angle difference (deg.) */
+#ifndef MAXANG
+#define MAXANG		20
+#endif
+#if MAXANG>0
+#define MAXDIFF2	( MAXANG*MAXANG * (PI*PI/180./180.))
+#endif
 
+#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 +34,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 +103,8 @@ newleaf()			/* allocate a leaf from our pile */
 }
 
 
-#define	LEAFSIZ		(3*sizeof(float)+sizeof(TMbright)+6*sizeof(BYTE))
+#define	LEAFSIZ		(3*sizeof(float)+sizeof(int4)+\
+			sizeof(TMbright)+6*sizeof(BYTE))
 
 int
 qtAllocLeaves(n)		/* allocate space for n leaves */
@@ -118,7 +130,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 = (int4 *)(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 +203,96 @@ int	pct;
 }
 
 
+#define	DCSCALE		11585.2		/* (1<<13)*sqrt(2) */
+#define FXNEG		01
+#define FYNEG		02
+#define FZNEG		04
+#define F1X		010
+#define F2Z		020
+#define F1SFT		5
+#define F2SFT		18
+#define FMASK		0x1fff
+
+static int4
+encodedir(dv)		/* encode a normalized direction vector */
+FVECT	dv;
+{
+	register int4	dc = 0;
+	int	cd[3], cm;
+	register int	i;
+
+	for (i = 0; i < 3; i++)
+		if (dv[i] < 0.) {
+			cd[i] = dv[i] * -DCSCALE;
+			dc |= 1<<i;
+		} else
+			cd[i] = dv[i] * DCSCALE;
+	if (cd[0] <= cd[1]) {
+		dc |= F1X | cd[0] << F1SFT;
+		cm = cd[1];
+	} else {
+		dc |= cd[1] << F1SFT;
+		cm = cd[0];
+	}
+	if (cd[2] <= cm)
+		dc |= F2Z | cd[2] << F2SFT;
+	else
+		dc |= cm << F2SFT;
+	return(dc);
+}
+
+
+static
+decodedir(dv, dc)	/* decode a normalized direction vector */
+register FVECT	dv;	/* returned */
+register int4	dc;
+{
+	double	d1, d2, der;
+
+	d1 = ((dc>>F1SFT & FMASK)+.5)/DCSCALE;
+	d2 = ((dc>>F2SFT & FMASK)+.5)/DCSCALE;
+	der = sqrt(1. - d1*d1 - d2*d2);
+	if (dc & F1X) {
+		dv[0] = d1;
+		if (dc & F2Z) { dv[1] = der; dv[2] = d2; }
+		else { dv[1] = d2; dv[2] = der; }
+	} else {
+		dv[1] = d1;
+		if (dc & F2Z) { dv[0] = der; dv[2] = d2; }
+		else { dv[0] = d2; dv[2] = der; }
+	}
+	if (dc & FXNEG) dv[0] = -dv[0];
+	if (dc & FYNEG) dv[1] = -dv[1];
+	if (dc & FZNEG) dv[2] = -dv[2];
+}
+
+
+static double
+dir2diff(dc1, dc2)		/* relative radians^2 between directions */
+int4	dc1, dc2;
+{
+	FVECT	v1, v2;
+
+	decodedir(v1, dc1);
+	decodedir(v2, dc2);
+
+	return(2. - 2.*DOT(v1,v2));
+}
+
+
+static double
+fdir2diff(dc1, v2)		/* relative radians^2 between directions */
+int4	dc1;
+register FVECT	v2;
+{
+	FVECT	v1;
+
+	decodedir(v1, dc1);
+
+	return(2. - 2.*DOT(v1,v2));
+}
+
+
 int
 qtFindLeaf(x, y)		/* find closest leaf to (x,y) */
 int	x, y;
@@ -234,16 +337,28 @@ int	li;
 	register RTREE	*tp = &qtrunk;
 	int	x0=0, y0=0, x1=odev.hres, y1=odev.vres;
 	int	lo = -1;
+	double	d2;
 	int	x, y, mx, my;
 	double	z;
-	FVECT	ip, wp;
+	FVECT	ip, wp, vd;
 	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(-1);			/* 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++)
+		vd[q] = (wp[q] - odev.v.vp[q])/ip[2];
+	d2 = fdir2diff(qtL.wd[li], vd);
+#ifdef MAXDIFF2
+	if (d2 > MAXDIFF2)
+		return(0);			/* leaf dir. too far off */
+#endif
 	x = ip[0] * odev.hres;
 	y = ip[1] * odev.vres;
 	z = ip[2];
@@ -266,17 +381,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] &&
+					fdir2diff(qtL.wd[lo], vd) < d2)
+				return(0);		/* old one better */
+			tp->k[q].li = li;		/* else new one is */
 			tp->flgs |= CHF(q);
 			break;
 		}
@@ -288,22 +405,25 @@ int	li;
 		my = ip[1] * odev.vres;
 		if (mx >= (x0 + x1) >> 1) q |= 01;
 		if (my >= (y0 + y1) >> 1) q |= 02;
+		tp->flgs = CH_ANY|LFF(q);	/* all new */
 		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);
+	qtL.wd[li] = encodedir(v);
 	tmCvColrs(&qtL.brt[li], qtL.chr[li], c, 1);
-	addleaf(li);
+	if (!addleaf(li))
+		ungetleaf(li);
 }