--- ray/src/hd/rhpict2.c	1999/03/08 14:09:11	3.3
+++ ray/src/hd/rhpict2.c	1999/03/09 15:08:22	3.8
@@ -10,6 +10,7 @@ static char SCCSid[] = "$SunId$ SGI";
 
 #include "holo.h"
 #include "view.h"
+#include "random.h"
 
 #ifndef DEPS
 #define DEPS		0.02	/* depth epsilon */
@@ -120,8 +121,12 @@ int	n;
 	int	d;
 	register int4	i, p;
 
-	if (n <= 0)
+	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--; ) {
@@ -139,7 +144,7 @@ int	n;
 
 
 int
-grow_samp(h, v, nl, n)		/* grow sample point appropriately */
+smooth_samp(h, v, nl, n)	/* grow sample point smoothly */
 int	h, v;
 register short	nl[NNEIGH][2];
 int	n;
@@ -147,7 +152,7 @@ int	n;
 	int	dis[NNEIGH], ndis;
 	COLOR	mykern[MAXRAD2];
 	int4	maxr2;
-	double	w, d;
+	double	d;
 	register int4	p, r2;
 	int	i, r, maxr, h2, v2;
 
@@ -192,34 +197,81 @@ int	n;
 			}
 			if (i >= 0) continue;	/* outside edge */
 			addcolor(mypixel[v2*hres+h2], mykern[r2]);
-			myweight[v2*hres+h2] += pixWeight[r2] *
-							myweight[v*hres+h];
+			myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
 		}
 	}
 	return(1);
 }
 
 
-pixFlush()			/* done with beams -- flush pixel values */
+int
+random_samp(h, v, nl, n, rf)	/* grow sample point noisily */
+int	h, v;
+register short	nl[NNEIGH][2];
+int	n;
+double	*rf;
 {
-	reset_flags();			/* set occupancy flags */
-	meet_neighbors(kill_occl);	/* eliminate occlusion errors */
-	reset_flags();			/* reset occupancy flags */
-	if (pixWeight[0] <= FTINY) {	/* initialize weighting function */
-		register int	i, j, r2;
-		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;
+	double	tv = *rf * (1.-1./NNEIGH);
+	int4	maxr2;
+	register int4	p, r2;
+	register int	i;
+	int	maxr, h2, v2;
+	COLOR	ctmp;
+
+	if (n <= 0)
+		return(1);
+	p = v*hres + h;				/* compute kernel radius */
+	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);
+						/* sample 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 */
+			if (frandom() < tv) {	/* use neighbor instead? */
+				i = random() % n;
+				r2 = nl[i][1]*hres + nl[i][0];
+				copycolor(ctmp, mypixel[r2]);
+				r2 = (h2-nl[i][0])*(h2-nl[i][0]) +
+						(v2-nl[i][1])*(v2-nl[i][1]);
+				if (r2 < MAXRAD2) {
+					scalecolor(ctmp, pixWeight[r2]);
+					addcolor(mypixel[v2*hres+h2], ctmp);
+					myweight[v2*hres+h2] += pixWeight[r2] *
+					      myweight[nl[i][1]*hres+nl[i][0]];
+					continue;
+				}
 			}
-		pixWeight[0] = 1.;
-		isqrttab[0] = 0;
+						/* use central sample */
+			copycolor(ctmp, mypixel[p]);
+			scalecolor(ctmp, pixWeight[r2]);
+			addcolor(mypixel[v2*hres+h2], ctmp);
+			myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
+		}
 	}
-	meet_neighbors(grow_samp);	/* grow valid samples over image */
+	return(1);
+}
+
+
+pixFinish(ransamp)		/* done with beams -- compute pixel values */
+double	ransamp;
+{
+	if (pixWeight[0] <= FTINY)
+		init_wfunc();		/* initialize weighting function */
+	reset_flags();			/* set occupancy flags */
+	meet_neighbors(kill_occl,NULL);	/* eliminate occlusion errors */
+	reset_flags();			/* reset occupancy flags */
+	if (ransamp >= 0.)		/* spread samples over image */
+		meet_neighbors(random_samp,&ransamp);
+	else
+		meet_neighbors(smooth_samp,NULL);
 	free((char *)pixFlags);		/* free pixel flags */
 	pixFlags = NULL;
 }
@@ -227,7 +279,7 @@ pixFlush()			/* done with beams -- flush pixel values 
 
 reset_flags()			/* allocate/set/reset occupancy flags */
 {
-	register int	p;
+	register int4	p;
 
 	if (pixFlags == NULL) {
 		pixFlags = (int4 *)calloc(FL4NELS(hres*vres), sizeof(int4));
@@ -240,6 +292,21 @@ reset_flags()			/* allocate/set/reset occupancy flags 
 }
 
 
+init_wfunc()			/* initialize weighting function */
+{
+	register int4	r2;
+	register double	d;
+
+	for (r2 = MAXRAD2; --r2; ) {
+		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];
@@ -281,8 +348,9 @@ register short	(*rnl)[NNEIGH];
 }
 
 
-meet_neighbors(nf)		/* run through samples and their neighbors */
+meet_neighbors(nf, dp)		/* run through samples and their neighbors */
 int	(*nf)();
+char	*dp;
 {
 	short	ln[NNEIGH][2];
 	int	h, v, n, v2;
@@ -306,7 +374,7 @@ int	(*nf)();
 			if (!CHK4(pixFlags, v*hres+h))
 				continue;	/* no one home */
 			n = findneigh(ln, h, v, rnl);
-			(*nf)(h, v, ln, n);	/* call on neighbors */
+			(*nf)(h, v, ln, n, dp);	/* call on neighbors */
 		}
 		if (++v >= vres)		/* reinitialize row list */
 			break;