--- ray/src/rt/ambcomp.c	2018/04/11 17:05:59	2.80
+++ ray/src/rt/ambcomp.c	2018/04/12 18:02:45	2.81
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: ambcomp.c,v 2.80 2018/04/11 17:05:59 greg Exp $";
+static const char	RCSid[] = "$Id: ambcomp.c,v 2.81 2018/04/12 18:02:45 greg Exp $";
 #endif
 /*
  * Routines to compute "ambient" values using Monte Carlo
@@ -156,7 +156,7 @@ resample:
 }
 
 
-/* Estimate errors based on ambient division differences */
+/* Estimate variance based on relative ambient division differences */
 static float *
 getambdiffs(AMBHEMI *hp)
 {
@@ -164,31 +164,34 @@ getambdiffs(AMBHEMI *hp)
 	float	*earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
 	float	*ep;
 	AMBSAMP	*ap;
-	double	b, d2;
+	double	b, b1, d2;
 	int	i, j;
 
 	if (earr == NULL)		/* out of memory? */
 		return(NULL);
-					/* compute squared neighbor diffs */
+					/* sum squared neighbor diffs */
 	for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
 	    for (j = 0; j < hp->ns; j++, ap++, ep++) {
 		b = bright(ap[0].v);
 		if (i) {		/* from above */
-			d2 = normf*(b - bright(ap[-hp->ns].v));
-			d2 *= d2;
+			b1 = bright(ap[-hp->ns].v);
+			d2 = (b - b1)/(b + b1);
+			d2 *= d2*normf;
 			ep[0] += d2;
 			ep[-hp->ns] += d2;
 		}
 		if (!j) continue;
 					/* from behind */
-		d2 = normf*(b - bright(ap[-1].v));
-		d2 *= d2;
+		b1 = bright(ap[-1].v);
+		d2 = (b - b1)/(b + b1);
+		d2 *= d2*normf;
 		ep[0] += d2;
 		ep[-1] += d2;
 		if (!i) continue;
 					/* diagonal */
-		d2 = normf*(b - bright(ap[-hp->ns-1].v));
-		d2 *= d2;
+		b1 = bright(ap[-hp->ns-1].v);
+		d2 = (b - b1)/(b + b1);
+		d2 *= d2*normf;
 		ep[0] += d2;
 		ep[-hp->ns-1] += d2;
 	    }
@@ -215,7 +218,6 @@ ambsupersamp(AMBHEMI *hp, int cnt)
 {
 	float	*earr = getambdiffs(hp);
 	double	e2rem = 0;
-	AMBSAMP	*ap;
 	float	*ep;
 	int	i, j, n, nss;
 
@@ -225,8 +227,8 @@ ambsupersamp(AMBHEMI *hp, int cnt)
 	for (ep = earr + hp->ns*hp->ns; ep > earr; )
 		e2rem += *--ep;
 	ep = earr;			/* perform super-sampling */
-	for (ap = hp->sa, i = 0; i < hp->ns; i++)
-	    for (j = 0; j < hp->ns; j++, ap++) {
+	for (i = 0; i < hp->ns; i++)
+	    for (j = 0; j < hp->ns; j++) {
 		if (e2rem <= FTINY)
 			goto done;	/* nothing left to do */
 		nss = *ep/e2rem*cnt + frandom();