--- ray/src/rt/ambcomp.c	2014/04/30 23:38:58	2.41
+++ ray/src/rt/ambcomp.c	2014/05/01 16:06:11	2.44
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: ambcomp.c,v 2.41 2014/04/30 23:38:58 greg Exp $";
+static const char	RCSid[] = "$Id: ambcomp.c,v 2.44 2014/05/01 16:06:11 greg Exp $";
 #endif
 /*
  * Routines to compute "ambient" values using Monte Carlo
@@ -22,20 +22,20 @@ static const char	RCSid[] = "$Id: ambcomp.c,v 2.41 201
 extern void		SDsquare2disk(double ds[2], double seedx, double seedy);
 
 typedef struct {
+	COLOR	v;		/* hemisphere sample value */
+	FVECT	p;		/* intersection point */
+} AMBSAMP;		/* sample value */
+
+typedef struct {
 	RAY	*rp;		/* originating ray sample */
 	FVECT	ux, uy;		/* tangent axis unit vectors */
 	int	ns;		/* number of samples per axis */
 	COLOR	acoef;		/* division contribution coefficient */
-	struct s_ambsamp {
-		COLOR	v;		/* hemisphere sample value */
-		FVECT	p;		/* intersection point */
-	} sa[1];		/* sample array (extends struct) */
+	AMBSAMP	sa[1];		/* sample array (extends struct) */
 }  AMBHEMI;		/* ambient sample hemisphere */
 
-typedef struct s_ambsamp	AMBSAMP;
+#define ambsam(h,i,j)	(h)->sa[(i)*(h)->ns + (j)]
 
-#define ambsamp(h,i,j)	(h)->sa[(i)*(h)->ns + (j)]
-
 typedef struct {
 	FVECT	r_i, r_i1, e_i, rcp, rI2_eJ2;
 	double	I1, I2;
@@ -88,9 +88,9 @@ inithemi(			/* initialize sampling hemisphere */
 }
 
 
-/* Prepare ambient division sample */
+/* Sample ambient division and apply weighting coefficient */
 static int
-prepambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
+getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
 {
 	int	hlist[3], ii;
 	double	spt[2], zd;
@@ -122,32 +122,32 @@ prepambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n
 				spt[1]*hp->uy[ii] +
 				zd*hp->rp->ron[ii];
 	checknorm(arp->rdir);
+	dimlist[ndims++] = i*hp->ns + j + 90171;
+	rayvalue(arp);			/* evaluate ray */
+	ndims--;			/* apply coefficient */
+	multcolor(arp->rcol, arp->rcoef);
 	return(1);
 }
 
 
 static AMBSAMP *
-ambsample(				/* sample an ambient direction */
+ambsample(				/* initial ambient division sample */
 	AMBHEMI	*hp,
 	int	i,
 	int	j
 )
 {
-	AMBSAMP	*ap = &ambsamp(hp,i,j);
+	AMBSAMP	*ap = &ambsam(hp,i,j);
 	RAY	ar;
 					/* generate hemispherical sample */
-	if (!prepambsamp(&ar, hp, i, j, 0))
+	if (!getambsamp(&ar, hp, i, j, 0))
 		goto badsample;
-	dimlist[ndims++] = i*hp->ns + j + 90171;
-	rayvalue(&ar);			/* evaluate ray */
-	ndims--;
 					/* limit vertex distance */
 	if (ar.rt > 10.0*thescene.cusize)
 		ar.rt = 10.0*thescene.cusize;
 	else if (ar.rt <= FTINY)	/* should never happen! */
 		goto badsample;
 	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
-	multcolor(ar.rcol, ar.rcoef);	/* apply coefficient */
 	copycolor(ap->v, ar.rcol);
 	return(ap);
 badsample:
@@ -163,23 +163,24 @@ getambdiffs(AMBHEMI *hp)
 {
 	float	*earr = calloc(hp->ns*hp->ns, sizeof(float));
 	float	*ep;
+	AMBSAMP	*ap;
 	double	b, d2;
 	int	i, j;
 
 	if (earr == NULL)		/* out of memory? */
 		return(NULL);
 					/* compute squared neighbor diffs */
-	for (ep = earr, i = 0; i < hp->ns; i++)
-	    for (j = 0; j < hp->ns; j++, ep++) {
-		b = bright(ambsamp(hp,i,j).v);
+	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 = b - bright(ambsamp(hp,i-1,j).v);
+			d2 = b - bright(ap[-hp->ns].v);
 			d2 *= d2;
 			ep[0] += d2;
 			ep[-hp->ns] += d2;
 		}
 		if (j) {		/* from behind */
-			d2 = b - bright(ambsamp(hp,i,j-1).v);
+			d2 = b - bright(ap[-1].v);
 			d2 *= d2;
 			ep[0] += d2;
 			ep[-1] += d2;
@@ -202,7 +203,7 @@ getambdiffs(AMBHEMI *hp)
 }
 
 
-/* Perform super-sampling on hemisphere */
+/* Perform super-sampling on hemisphere (introduces bias) */
 static void
 ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
 {
@@ -225,14 +226,10 @@ ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
 		int	nss = *ep/e2sum*cnt + frandom();
 		setcolor(asum, 0., 0., 0.);
 		for (n = 1; n <= nss; n++) {
-			if (!prepambsamp(&ar, hp, i, j, n)) {
+			if (!getambsamp(&ar, hp, i, j, n)) {
 				nss = n-1;
 				break;
 			}
-			dimlist[ndims++] = i*hp->ns + j + 90171;
-			rayvalue(&ar);	/* evaluate super-sample */
-			ndims--;
-			multcolor(ar.rcol, ar.rcoef);
 			addcolor(asum, ar.rcol);
 		}
 		if (nss) {		/* update returned ambient value */
@@ -495,8 +492,8 @@ ambHessian(				/* anisotropic radii & pos. gradient */
 	}
 					/* compute first row of edges */
 	for (j = 0; j < hp->ns-1; j++) {
-		comp_fftri(&fftr, ambsamp(hp,0,j).p,
-				ambsamp(hp,0,j+1).p, hp->rp->rop);
+		comp_fftri(&fftr, ambsam(hp,0,j).p,
+				ambsam(hp,0,j+1).p, hp->rp->rop);
 		if (hessrow != NULL)
 			comp_hessian(hessrow[j], &fftr, hp->rp->ron);
 		if (gradrow != NULL)
@@ -506,8 +503,8 @@ ambHessian(				/* anisotropic radii & pos. gradient */
 	for (i = 0; i < hp->ns-1; i++) {
 	    FVECT	hesscol[3];	/* compute first vertical edge */
 	    FVECT	gradcol;
-	    comp_fftri(&fftr, ambsamp(hp,i,0).p,
-			ambsamp(hp,i+1,0).p, hp->rp->rop);
+	    comp_fftri(&fftr, ambsam(hp,i,0).p,
+			ambsam(hp,i+1,0).p, hp->rp->rop);
 	    if (hessrow != NULL)
 		comp_hessian(hesscol, &fftr, hp->rp->ron);
 	    if (gradrow != NULL)
@@ -516,11 +513,11 @@ ambHessian(				/* anisotropic radii & pos. gradient */
 		FVECT	hessdia[3];	/* compute triangle contributions */
 		FVECT	graddia;
 		COLORV	backg;
-		backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
-					&ambsamp(hp,i+1,j), hp->rp->rop);
+		backg = back_ambval(&ambsam(hp,i,j), &ambsam(hp,i,j+1),
+					&ambsam(hp,i+1,j), hp->rp->rop);
 					/* diagonal (inner) edge */
-		comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
-				ambsamp(hp,i+1,j).p, hp->rp->rop);
+		comp_fftri(&fftr, ambsam(hp,i,j+1).p,
+				ambsam(hp,i+1,j).p, hp->rp->rop);
 		if (hessrow != NULL) {
 		    comp_hessian(hessdia, &fftr, hp->rp->ron);
 		    rev_hessian(hesscol);
@@ -532,16 +529,16 @@ ambHessian(				/* anisotropic radii & pos. gradient */
 		    add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
 		}
 					/* initialize edge in next row */
-		comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
-				ambsamp(hp,i+1,j).p, hp->rp->rop);
+		comp_fftri(&fftr, ambsam(hp,i+1,j+1).p,
+				ambsam(hp,i+1,j).p, hp->rp->rop);
 		if (hessrow != NULL)
 		    comp_hessian(hessrow[j], &fftr, hp->rp->ron);
 		if (gradrow != NULL)
 		    comp_gradient(gradrow[j], &fftr, hp->rp->ron);
 					/* new column edge & paired triangle */
-		backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
-					&ambsamp(hp,i+1,j), hp->rp->rop);
-		comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
+		backg = back_ambval(&ambsam(hp,i,j+1), &ambsam(hp,i+1,j+1),
+					&ambsam(hp,i+1,j), hp->rp->rop);
+		comp_fftri(&fftr, ambsam(hp,i,j+1).p, ambsam(hp,i+1,j+1).p,
 				hp->rp->rop);
 		if (hessrow != NULL) {
 		    comp_hessian(hesscol, &fftr, hp->rp->ron);