--- ray/src/hd/rhpict2.c	1999/03/04 10:30:04	3.1
+++ ray/src/hd/rhpict2.c	1999/03/09 11:41:29	3.5
@@ -11,13 +11,49 @@ static char SCCSid[] = "$SunId$ SGI";
 #include "holo.h"
 #include "view.h"
 
+#ifndef DEPS
+#define DEPS		0.02	/* depth epsilon */
+#endif
+#ifndef PEPS
+#define PEPS		0.04	/* pixel value epsilon */
+#endif
+#ifndef MAXRAD
+#define MAXRAD		64	/* maximum kernel radius */
+#endif
+#ifndef NNEIGH
+#define NNEIGH		7	/* find this many neighbors */
+#endif
+
+#define NINF		16382
+
+#define MAXRAD2		(MAXRAD*MAXRAD+1)
+
+#define G0NORM		0.286	/* ground zero normalization (1/x integral) */
+
+#ifndef FL4OP
+#define FL4OP(f,i,op)	((f)[(i)>>5] op (1L<<((i)&0x1f)))
+#define CHK4(f,i)	FL4OP(f,i,&)
+#define SET4(f,i)	FL4OP(f,i,|=)
+#define CLR4(f,i)	FL4OP(f,i,&=~)
+#define TGL4(f,i)	FL4OP(f,i,^=)
+#define FL4NELS(n)	(((n)+0x1f)>>5)
+#define CLR4ALL(f,n)	bzero((char *)(f),FL4NELS(n)*sizeof(int4))
+#endif
+
+static int4	*pixFlags;	/* pixel occupancy flags */
+static float	pixWeight[MAXRAD2];	/* pixel weighting function */
+static short	isqrttab[MAXRAD2];	/* integer square root table */
+
+#define isqrt(i2)	((int)isqrttab[(int)(i2)])
+
 extern VIEW	myview;		/* current output view */
 extern COLOR	*mypixel;	/* pixels being rendered */
 extern float	*myweight;	/* weights (used to compute final pixels) */
+extern float	*mydepth;	/* depth values (visibility culling) */
 extern int	hres, vres;	/* current horizontal and vertical res. */
 
 
-render_beam(bp, hb)		/* render a particular beam */
+pixBeam(bp, hb)			/* render a particular beam */
 BEAM	*bp;
 register HDBEAMI	*hb;
 {
@@ -26,7 +62,8 @@ register HDBEAMI	*hb;
 	FVECT	rorg, rdir, wp, ip;
 	double	d, prox;
 	COLOR	col;
-	int	n, p;
+	int	n;
+	register int4	p;
 
 	if (!hdbcoord(gc, hb->h, hb->b))
 		error(CONSISTENCY, "bad beam in render_beam");
@@ -38,12 +75,13 @@ register HDBEAMI	*hb;
 			VSUM(wp, rorg, rdir, d);
 			VSUB(ip, wp, myview.vp);
 			d = DOT(ip,rdir);
-			prox = 1. - d*d/DOT(ip,ip);	/* sin(diff_angle)^2 */
+			prox = d*d/DOT(ip,ip);	/* cos(diff_angle)^32 */
+			prox *= prox; prox *= prox; prox *= prox; prox *= prox;
 		} else {
 			if (myview.type == VT_PAR || myview.vaft > FTINY)
 				continue;		/* inf. off view */
 			VSUM(wp, myview.vp, rdir, FHUGE);
-			prox = 0.;
+			prox = 1.;
 		}
 		viewloc(ip, &myview, wp);	/* frustum clipping */
 		if (ip[2] < 0.)
@@ -53,11 +91,245 @@ register HDBEAMI	*hb;
 		if (ip[1] < 0. || ip[1] >= 1.)
 			continue;
 		if (myview.vaft > FTINY && ip[2] > myview.vaft - myview.vfore)
-			continue;
-						/* we're in! */
+			continue;		/* not exact for VT_PER */
 		p = (int)(ip[1]*vres)*hres + (int)(ip[0]*hres);
+		if (mydepth[p] > FTINY) {	/* check depth */
+			if (ip[2] > mydepth[p]*(1.+DEPS))
+				continue;
+			if (ip[2] < mydepth[p]*(1.-DEPS)) {
+				setcolor(mypixel[p], 0., 0., 0.);
+				myweight[p] = 0.;
+			}
+		}
 		colr_color(col, rv->v);
+		scalecolor(col, prox);
 		addcolor(mypixel[p], col);
-		myweight[p] += 1.0;
+		myweight[p] += prox;
+		mydepth[p] = ip[2];
 	}
+}
+
+
+int
+kill_occl(h, v, nl, n)		/* check for occlusion errors */
+int	h, v;
+register short	nl[NNEIGH][2];
+int	n;
+{
+	short	forequad[2][2];
+	int	d;
+	register int4	i, p;
+
+	if (n <= 0) {
+#ifdef DEBUG
+		error(WARNING, "neighborless sample in kill_occl");
+#endif
+		return(1);
+	}
+	p = v*hres + h;
+	forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
+	for (i = n; i--; ) {
+		d = (h-nl[i][0])*(h-nl[i][0]) + (v-nl[i][1])*(v-nl[i][1]);
+		d = isqrt(d);
+		if (mydepth[nl[i][1]*hres+nl[i][0]]*(1.+DEPS*d) < mydepth[p])
+			forequad[nl[i][0]<h][nl[i][1]<v] = 1;
+	}
+	if (forequad[0][0]+forequad[0][1]+forequad[1][0]+forequad[1][1] > 2) {
+		setcolor(mypixel[p], 0., 0., 0.);
+		myweight[p] = 0.;	/* occupancy reset afterwards */
+	}
+	return(1);
+}
+
+
+int
+grow_samp(h, v, nl, n)		/* grow sample point appropriately */
+int	h, v;
+register short	nl[NNEIGH][2];
+int	n;
+{
+	int	dis[NNEIGH], ndis;
+	COLOR	mykern[MAXRAD2];
+	int4	maxr2;
+	double	d;
+	register int4	p, r2;
+	int	i, r, maxr, h2, v2;
+
+	if (n <= 0)
+		return(1);
+	p = v*hres + h;				/* build kernel values */
+	maxr2 = (h-nl[n-1][0])*(h-nl[n-1][0]) + (v-nl[n-1][1])*(v-nl[n-1][1]);
+	DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
+	maxr = isqrt(maxr2);
+	for (v2 = 1; v2 <= maxr; v2++)
+		for (h2 = 0; h2 <= v2; h2++) {
+			r2 = h2*h2 + v2*v2;
+			if (r2 > maxr2) break;
+			copycolor(mykern[r2], mypixel[p]);
+			scalecolor(mykern[r2], pixWeight[r2]);
+		}
+	ndis = 0;				/* find discontinuities */
+	for (i = n; i--; ) {
+		r2 = (h-nl[i][0])*(h-nl[i][0]) + (v-nl[i][1])*(v-nl[i][1]);
+		r = isqrt(r2);
+		d = mydepth[nl[i][1]*hres+nl[i][0]] / mydepth[p];
+		d = d>=1. ? d-1. : 1.-d;
+		if (d > r*DEPS || bigdiff(mypixel[p],
+				mypixel[nl[i][1]*hres+nl[i][0]], r*PEPS))
+			dis[ndis++] = i;
+	}
+						/* stamp out that kernel */
+	for (v2 = v-maxr; v2 <= v+maxr; v2++) {
+		if (v2 < 0) v2 = 0;
+		else if (v2 >= vres) break;
+		for (h2 = h-maxr; h2 <= h+maxr; h2++) {
+			if (h2 < 0) h2 = 0;
+			else if (h2 >= hres) break;
+			r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
+			if (r2 > maxr2) continue;
+			if (CHK4(pixFlags, v2*hres+h2))
+				continue;	/* occupied */
+			for (i = ndis; i--; ) {
+				r = (h2-nl[dis[i]][0])*(h2-nl[dis[i]][0]) +
+					(v2-nl[dis[i]][1])*(v2-nl[dis[i]][1]);
+				if (r < r2) break;
+			}
+			if (i >= 0) continue;	/* outside edge */
+			addcolor(mypixel[v2*hres+h2], mykern[r2]);
+			myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
+		}
+	}
+	return(1);
+}
+
+
+pixFlush()			/* done with beams -- flush pixel values */
+{
+	if (pixWeight[0] <= FTINY)
+		init_wfunc();		/* initialize weighting function */
+	reset_flags();			/* set occupancy flags */
+	meet_neighbors(kill_occl);	/* eliminate occlusion errors */
+	reset_flags();			/* reset occupancy flags */
+	meet_neighbors(grow_samp);	/* grow valid samples over image */
+	free((char *)pixFlags);		/* free pixel flags */
+	pixFlags = NULL;
+}
+
+
+reset_flags()			/* allocate/set/reset occupancy flags */
+{
+	register int4	p;
+
+	if (pixFlags == NULL) {
+		pixFlags = (int4 *)calloc(FL4NELS(hres*vres), sizeof(int4));
+		CHECK(pixFlags==NULL, SYSTEM, "out of memory in reset_flags");
+	} else
+		CLR4ALL(pixFlags, hres*vres);
+	for (p = hres*vres; p--; )
+		if (myweight[p] > FTINY)
+			SET4(pixFlags, p);
+}
+
+
+init_wfunc()			/* initialize weighting function */
+{
+	register int	i, j;
+	register int4	r2;
+	register double	d;
+
+	for (i = 1; i <= MAXRAD; i++)
+		for (j = 0; j <= i; j++) {
+			r2 = i*i + j*j;
+			if (r2 >= MAXRAD2) break;
+			d = sqrt((double)r2);
+			pixWeight[r2] = G0NORM/d;
+			isqrttab[r2] = d + 0.99;
+		}
+	pixWeight[0] = 1.;
+	isqrttab[0] = 0;
+}
+
+
+int
+findneigh(nl, h, v, rnl)	/* find NNEIGH neighbors for pixel */
+short	nl[NNEIGH][2];
+int	h, v;
+register short	(*rnl)[NNEIGH];
+{
+	int	nn = 0;
+	int4	d, nd[NNEIGH];
+	int	n, hoff;
+	register int	h2, n2;
+
+	nd[NNEIGH-1] = MAXRAD2;
+	for (hoff = 1; hoff < hres; hoff = (hoff<0) - hoff) {
+		h2 = h + hoff;
+		if (h2 < 0 | h2 >= hres)
+			continue;
+		if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
+			break;
+		for (n = 0; n < NNEIGH && rnl[h2][n] < NINF; n++) {
+			d = (h2-h)*(h2-h) + (v-rnl[h2][n])*(v-rnl[h2][n]);
+			if (d >= nd[NNEIGH-1])
+				continue;
+			if (nn < NNEIGH)	/* insert neighbor */
+				nn++;
+			for (n2 = nn; n2--; ) 	{
+				if (!n2 || d >= nd[n2-1]) {
+					nd[n2] = d;
+					nl[n2][0] = h2;
+					nl[n2][1] = rnl[h2][n];
+					break;
+				}
+				nd[n2] = nd[n2-1];
+				nl[n2][0] = nl[n2-1][0];
+				nl[n2][1] = nl[n2-1][1];
+			}
+		}
+	}
+	return(nn);
+}
+
+
+meet_neighbors(nf)		/* run through samples and their neighbors */
+int	(*nf)();
+{
+	short	ln[NNEIGH][2];
+	int	h, v, n, v2;
+	register short	(*rnl)[NNEIGH];
+					/* initialize bottom row list */
+	rnl = (short (*)[NNEIGH])malloc(NNEIGH*sizeof(short)*hres);
+	CHECK(rnl==NULL, SYSTEM, "out of memory in meet_neighbors");
+	for (h = 0; h < hres; h++) {
+		for (n = v = 0; v < vres; v++)
+			if (CHK4(pixFlags, v*hres+h)) {
+				rnl[h][n++] = v;
+				if (n >= NNEIGH)
+					break;
+			}
+		while (n < NNEIGH)
+			rnl[h][n++] = NINF;
+	}
+	v = 0;				/* do each row */
+	for ( ; ; ) {
+		for (h = 0; h < hres; h++) {
+			if (!CHK4(pixFlags, v*hres+h))
+				continue;	/* no one home */
+			n = findneigh(ln, h, v, rnl);
+			(*nf)(h, v, ln, n);	/* call on neighbors */
+		}
+		if (++v >= vres)		/* reinitialize row list */
+			break;
+		for (h = 0; h < hres; h++)
+			for (v2 = rnl[h][NNEIGH-1]+1; v2 < vres; v2++) {
+				if (v2 - v > v - rnl[h][0])
+					break;		/* not close enough */
+				if (CHK4(pixFlags, v2*hres+h)) {
+					for (n = 0; n < NNEIGH-1; n++)
+						rnl[h][n] = rnl[h][n+1];
+					rnl[h][NNEIGH-1] = v2;
+				}
+			}
+	}
+	free((char *)rnl);		/* free row list */
 }