--- ray/src/hd/rholo2.c	1998/11/24 17:38:09	3.17
+++ ray/src/hd/rholo2.c	1998/12/03 15:21:05	3.20
@@ -15,72 +15,178 @@ static char SCCSid[] = "$SunId$ SGI";
 
 VIEWPOINT	myeye;		/* target view position */
 
+struct gclim {
+	HOLO	*hp;			/* holodeck pointer */
+	GCOORD	gc;			/* grid cell */
+	FVECT	egp;			/* eye grid point */
+	double	erg2;			/* mean square eye grid range */
+	double	gmin[2], gmax[2];	/* grid coordinate limits */
+};				/* a grid coordinate range */
 
+
+static
+initeyelim(gcl, hp, gc)		/* initialize grid coordinate limits */
+register struct gclim	*gcl;
+register HOLO	*hp;
+GCOORD	*gc;
+{
+	register FLOAT	*v;
+	register int	i;
+
+	if (hp != NULL) {
+		hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
+		gcl->erg2 = 0;
+		for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
+			gcl->erg2 += DOT(v,v);
+		gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
+	}
+	if (gc != NULL)
+		copystruct(&gcl->gc, gc);
+	gcl->gmin[0] = gcl->gmin[1] = FHUGE;
+	gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
+}
+
+
+static
+groweyelim(gcl, gc, r0, r1)	/* grow grid limits about eye point */
+register struct gclim	*gcl;
+GCOORD	*gc;
+double	r0, r1;
+{
+	FVECT	gp, ab;
+	double	vlen, plen, dv0, dv1;
+	double	rd2, dwall, gpos;
+	int	eyeout;
+	register int	i, g0, g1;
+
+	i = gc->w>>1;
+	if (gc->w&1)
+		eyeout = (gp[i] = gcl->hp->grid[i]) < gcl->egp[i];
+	else
+		eyeout = (gp[i] = 0) > gcl->egp[i];
+	gp[hdwg0[gc->w]] = gc->i[0] + r0;
+	gp[hdwg1[gc->w]] = gc->i[1] + r1;
+	VSUB(ab, gcl->egp, gp);
+	rd2 = DOT(ab,ab);
+	if (rd2 <= gcl->erg2) {
+		gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
+		gcl->gmax[0] = gcl->gmax[1] = FHUGE;
+		return;
+	}
+	rd2 = gcl->erg2 / rd2;
+	vlen = 1. - rd2;
+	plen = sqrt(rd2 * vlen);
+	g0 = gcl->gc.w>>1;
+	dwall = (gcl->gc.w&1 ? gcl->hp->grid[g0] : 0) - gp[g0];
+	for (i = 0; i < 4; i++) {
+		if (i == 2)
+			plen = -plen;
+		g1 = (g0+(i&1)+1)%3;
+		dv0 = vlen*ab[g0] + plen*ab[g1];
+		dv1 = vlen*ab[g1] - plen*ab[g0];
+		if ((dv0 < 0 ^ dwall < 0 ^ eyeout) ||
+				(dv0 <= FTINY && dv0 >= -FTINY)) {
+			if (eyeout)
+				dv1 = -dv1;
+			if (dv1 > FTINY)
+				gcl->gmax[i&1] = FHUGE;
+			else if (dv1 < -FTINY)
+				gcl->gmin[i&1] = -FHUGE;
+		} else {
+			gpos = gp[g1] + dv1*dwall/dv0;
+			if (gpos < gcl->gmin[i&1])
+				gcl->gmin[i&1] = gpos;
+			if (gpos > gcl->gmax[i&1])
+				gcl->gmax[i&1] = gpos;
+		}
+	}
+}
+
+
+static int
+clipeyelim(rrng, gcl)		/* clip eye limits to grid cell */
+register short	rrng[2][2];
+register struct gclim	*gcl;
+{
+	int	incell = 1;
+	register int	i;
+
+	for (i = 0; i < 2; i++) {
+		if (gcl->gmin[i] < gcl->gc.i[i])
+			gcl->gmin[i] = gcl->gc.i[i];
+		if (gcl->gmax[i] > gcl->gc.i[i]+1)
+			gcl->gmax[i] = gcl->gc.i[i]+1;
+		if (gcl->gmax[i] > gcl->gmin[i]) {
+			rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
+					(1.-FTINY);
+			rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
+					(1.-FTINY) - rrng[i][0];
+		} else
+			rrng[i][0] = rrng[i][1] = 0;
+		incell &= rrng[i][1] > 0;
+	}
+	return(incell);
+}
+
+
 packrays(rod, p)		/* pack ray origins and directions */
 register float	*rod;
 register PACKET	*p;
 {
-	short	packord[RPACKSIZ];
-	float	packdc2[RPACKSIZ];
-	int	iterleft = 3*p->nr;
-	BYTE	rpos[2][2];
-	FVECT	ro, rd, rp1;
+	int	nretries = p->nr + 2;
+	struct gclim	eyelim;
+	short	rrng0[2][2], rrng1[2][2];
+	int	useyelim;
 	GCOORD	gc[2];
-	double	d, dc2, md2, td2;
-	int	i;
-	register int	ii;
+	FVECT	ro, rd;
+	double	d;
+	register int	i;
 
 	if (!hdbcoord(gc, hdlist[p->hd], p->bi))
 		error(CONSISTENCY, "bad beam index in packrays");
-	td2 = (myeye.rng+FTINY)*(myeye.rng+FTINY);
-	for (i = 0, md2 = 0.; i < p->nr || md2 > td2; ) {
-		rpos[0][0] = frandom() * 256.;
-		rpos[0][1] = frandom() * 256.;
-		rpos[1][0] = frandom() * 256.;
-		rpos[1][1] = frandom() * 256.;
-		d = hdray(ro, rd, hdlist[p->hd], gc, rpos);
-		if (myeye.rng > FTINY) {		/* check eyepoint */
-			register int	nexti;
-
-			VSUM(rp1, ro, rd, d);
-			dc2 = dist2line(myeye.vpt, ro, rp1);
-			dc2 /= (double)(p->nr*p->nr);
-			if (i == p->nr) {		/* packet full */
-				nexti = packord[i-1];
-				if (!iterleft--)
-					break;		/* tried enough! */
-				if (dc2 >= packdc2[nexti])
-					continue;	/* worse than worst */
-				md2 -= packdc2[nexti];
-			} else
-				nexti = i++;
-			md2 += packdc2[nexti] = dc2;	/* new distance */
-			for (ii = i; --ii; ) {		/* insertion sort */
-				if (dc2 > packdc2[packord[ii-1]])
-					break;
-				packord[ii] = packord[ii-1];
+	if ((useyelim = myeye.rng > FTINY)) {
+		initeyelim(&eyelim, hdlist[p->hd], gc);
+		groweyelim(&eyelim, gc+1, 0., 0.);
+		groweyelim(&eyelim, gc+1, 1., 1.);
+		useyelim &= clipeyelim(rrng0, &eyelim);
+	}
+	for (i = 0; i < p->nr; i++) {
+	retry:
+		if (useyelim) {
+			initeyelim(&eyelim, NULL, gc+1);
+			p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
+						+ rrng0[0][0];
+			p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
+						+ rrng0[1][0];
+			groweyelim(&eyelim, gc,
+					(1./256.)*(p->ra[i].r[0][0]+.5),
+					(1./256.)*(p->ra[i].r[0][1]+.5));
+			if (!clipeyelim(rrng1, &eyelim)) {
+				useyelim &= nretries-- > 0;
+				goto retry;
 			}
-			packord[ii] = nexti;
-			ii = nexti;			/* put it here */
-		} else
-			ii = i++;
+			p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
+						+ rrng1[0][0];
+			p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
+						+ rrng1[1][0];
+		} else {
+			p->ra[i].r[0][0] = frandom() * 256.;
+			p->ra[i].r[0][1] = frandom() * 256.;
+			p->ra[i].r[1][0] = frandom() * 256.;
+			p->ra[i].r[1][1] = frandom() * 256.;
+		}
+		d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
 		if (p->offset != NULL) {
 			if (!vdef(OBSTRUCTIONS))
 				d *= frandom();		/* random offset */
 			VSUM(ro, ro, rd, d);		/* advance ray */
-			p->offset[ii] = d;
+			p->offset[i] = d;
 		}
-		p->ra[ii].r[0][0] = rpos[0][0];
-		p->ra[ii].r[0][1] = rpos[0][1];
-		p->ra[ii].r[1][0] = rpos[1][0];
-		p->ra[ii].r[1][1] = rpos[1][1];
-		VCOPY(rod+6*ii, ro);
-		VCOPY(rod+6*ii+3, rd);
+		VCOPY(rod, ro);
+		rod += 3;
+		VCOPY(rod, rd);
+		rod += 3;
 	}
-#ifdef DEBUG
-	fprintf(stderr, "%f mean distance for target %f (%d iterations left)\n",
-			sqrt(md2), myeye.rng, iterleft);
-#endif
 }