--- ray/src/rt/ambcomp.c	2014/04/19 19:20:47	2.28
+++ ray/src/rt/ambcomp.c	2014/04/24 06:03:15	2.31
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: ambcomp.c,v 2.28 2014/04/19 19:20:47 greg Exp $";
+static const char	RCSid[] = "$Id: ambcomp.c,v 2.31 2014/04/24 06:03:15 greg Exp $";
 #endif
 /*
  * Routines to compute "ambient" values using Monte Carlo
@@ -28,15 +28,15 @@ typedef struct {
 	COLOR	acoef;		/* division contribution coefficient */
 	struct s_ambsamp {
 		COLOR	v;		/* hemisphere sample value */
-		float	p[3];		/* intersection point */
+		FVECT	p;		/* intersection point */
 	} sa[1];		/* sample array (extends struct) */
 }  AMBHEMI;		/* ambient sample hemisphere */
 
 #define ambsamp(h,i,j)	(h)->sa[(i)*(h)->ns + (j)]
 
 typedef struct {
-	FVECT	r_i, r_i1, e_i;
-	double	nf, I1, I2, J2;
+	FVECT	r_i, r_i1, e_i, rI2_eJ2;
+	double	nf, I1, I2;
 } FFTRI;		/* vectors and coefficients for Hessian calculation */
 
 
@@ -96,7 +96,6 @@ ambsample(				/* sample an ambient direction */
 {
 	struct s_ambsamp	*ap = &ambsamp(hp,i,j);
 	RAY			ar;
-	int			hlist[3];
 	double			spt[2], zd;
 	int			ii;
 					/* ambient coefficient for weight */
@@ -104,11 +103,8 @@ ambsample(				/* sample an ambient direction */
 		setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
 	else
 		copycolor(ar.rcoef, hp->acoef);
-	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) {
-		setcolor(ap->v, 0., 0., 0.);
-		VCOPY(ap->p, hp->rp->rop);
-		return(NULL);		/* no sample taken */
-	}
+	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
+		goto badsample;
 	if (ambacc > FTINY) {
 		multcolor(ar.rcoef, hp->acoef);
 		scalecolor(ar.rcoef, 1./AVGREFL);
@@ -125,37 +121,45 @@ ambsample(				/* sample an ambient direction */
 	dimlist[ndims++] = i*hp->ns + j + 90171;
 	rayvalue(&ar);			/* evaluate ray */
 	ndims--;
-	multcolor(ar.rcol, ar.rcoef);	/* apply coefficient */
-	copycolor(ap->v, ar.rcol);
 	if (ar.rt > 20.0*maxarad)	/* limit vertex distance */
 		ar.rt = 20.0*maxarad;
+	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:
+	setcolor(ap->v, 0., 0., 0.);
+	VCOPY(ap->p, hp->rp->rop);
+	return(NULL);
 }
 
 
 /* Compute vectors and coefficients for Hessian/gradient calcs */
 static void
-comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop)
+comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop)
 {
-	FVECT	v1;
-	double	dot_e, dot_er, dot_r, dot_r1;
+	FVECT	vcp;
+	double	dot_e, dot_er, dot_r, dot_r1, J2;
+	int	i;
 
 	VSUB(ftp->r_i, ap0, rop);
 	VSUB(ftp->r_i1, ap1, rop);
 	VSUB(ftp->e_i, ap1, ap0);
-	VCROSS(v1, ftp->e_i, ftp->r_i);
-	ftp->nf = 1.0/DOT(v1,v1);
-	VCROSS(v1, ftp->r_i, ftp->r_i1);
-	ftp->I1 = sqrt(DOT(v1,v1)*ftp->nf);
+	VCROSS(vcp, ftp->e_i, ftp->r_i);
+	ftp->nf = 1.0/DOT(vcp,vcp);
 	dot_e = DOT(ftp->e_i,ftp->e_i);
 	dot_er = DOT(ftp->e_i, ftp->r_i);
 	dot_r = DOT(ftp->r_i,ftp->r_i);
 	dot_r1 = DOT(ftp->r_i1,ftp->r_i1);
+	ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) / sqrt(dot_r*dot_r1) ) *
+			sqrt( ftp->nf );
 	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r +
 			dot_e*ftp->I1 )*0.5*ftp->nf;
-	ftp->J2 =  0.25*ftp->nf*( 1.0/dot_r - 1.0/dot_r1 ) -
-			dot_er/dot_e*ftp->I2;
+	J2 =  0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - dot_er/dot_e*ftp->I2;
+	for (i = 3; i--; )
+		ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i];
 }
 
 
@@ -176,7 +180,7 @@ compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
 static void
 comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
 {
-	FVECT	v1, v2;
+	FVECT	vcp;
 	FVECT	m1[3], m2[3], m3[3], m4[3];
 	double	d1, d2, d3, d4;
 	double	I3, J3, K3;
@@ -187,19 +191,17 @@ comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
 	d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
 	d4 = DOT(ftp->e_i, ftp->r_i);
 	I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 +
-				3.0*d3*ftp->I2 );
+				3.0/d3*ftp->I2 );
 	J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3;
 	K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3);
 					/* intermediate matrices */
-	VCROSS(v1, nrm, ftp->e_i);
-	for (j = 3; j--; )
-		v2[j] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j];
-	compose_matrix(m1, v1, v2);
+	VCROSS(vcp, nrm, ftp->e_i);
+	compose_matrix(m1, vcp, ftp->rI2_eJ2);
 	compose_matrix(m2, ftp->r_i, ftp->r_i);
 	compose_matrix(m3, ftp->e_i, ftp->e_i);
 	compose_matrix(m4, ftp->r_i, ftp->e_i);
-	VCROSS(v1, ftp->r_i, ftp->e_i);
-	d1 = DOT(nrm, v1);
+	VCROSS(vcp, ftp->r_i, ftp->e_i);
+	d1 = DOT(nrm, vcp);
 	d2 = -d1*ftp->I2;
 	d1 *= 2.0;
 	for (i = 3; i--; )		/* final matrix sum */
@@ -251,8 +253,7 @@ comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
 	f1 = 2.0*DOT(nrm, vcp);
 	VCROSS(vcp, nrm, ftp->e_i);
 	for (i = 3; i--; )
-		grad[i] = (0.5/PI)*( ftp->I1*vcp[i] +
-			    f1*(ftp->I2*ftp->r_i[i] + ftp->J2*ftp->e_i[i]) );
+		grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] );
 }
 
 
@@ -288,17 +289,17 @@ back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *a
 
 	VSUB(vec, ap1->p, orig);
 	d2best = DOT(vec,vec);
-	vback = ap1->v[CIEY];
+	vback = colval(ap1->v,CIEY);
 	VSUB(vec, ap2->p, orig);
 	d2 = DOT(vec,vec);
 	if (d2 > d2best) {
 		d2best = d2;
-		vback = ap2->v[CIEY];
+		vback = colval(ap2->v,CIEY);
 	}
 	VSUB(vec, ap3->p, orig);
 	d2 = DOT(vec,vec);
 	if (d2 > d2best)
-		return(ap3->v[CIEY]);
+		return(colval(ap3->v,CIEY));
 	return(vback);
 }
 
@@ -328,12 +329,12 @@ eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
 		error(INTERNAL, "bad eigenvalue calculation");
 
 	if (evalue[0] > evalue[1]) {
-		ra[0] = 1.0/sqrt(sqrt(evalue[0]));
-		ra[1] = 1.0/sqrt(sqrt(evalue[1]));
+		ra[0] = sqrt(sqrt(4.0/evalue[0]));
+		ra[1] = sqrt(sqrt(4.0/evalue[1]));
 		slope1 = evalue[1];
 	} else {
-		ra[0] = 1.0/sqrt(sqrt(evalue[1]));
-		ra[1] = 1.0/sqrt(sqrt(evalue[0]));
+		ra[0] = sqrt(sqrt(4.0/evalue[1]));
+		ra[1] = sqrt(sqrt(4.0/evalue[0]));
 		slope1 = evalue[0];
 	}
 					/* compute unit eigenvectors */
@@ -454,9 +455,9 @@ ambHessian(				/* anisotropic radii & pos. gradient */
 	
 	if (ra != NULL)			/* extract eigenvectors & radii */
 		eigenvectors(uv, ra, hessian);
-	if (pg != NULL) {		/* project position gradient */
-		pg[0] = DOT(gradient, uv[0]);
-		pg[1] = DOT(gradient, uv[1]);
+	if (pg != NULL) {		/* tangential position gradient/PI */
+		pg[0] = DOT(gradient, uv[0]) / PI;
+		pg[1] = DOT(gradient, uv[1]) / PI;
 	}
 }
 
@@ -466,20 +467,23 @@ static void
 ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
 {
 	struct s_ambsamp	*ap;
+	double			dgsum[2];
 	int			n;
 	FVECT			vd;
 	double			gfact;
 
-	dg[0] = dg[1] = 0;
+	dgsum[0] = dgsum[1] = 0.0;	/* sum values times -tan(theta) */
 	for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
 					/* use vector for azimuth + 90deg */
 		VSUB(vd, ap->p, hp->rp->rop);
-					/* brightness with tangent factor */
-		gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd);
-					/* sine = proj_radius/vd_length */
-		dg[0] -= DOT(uv[1], vd) * gfact;
-		dg[1] += DOT(uv[0], vd) * gfact;
+					/* brightness over cosine factor */
+		gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
+					/* -sine = -proj_radius/vd_length */
+		dgsum[0] += DOT(uv[1], vd) * gfact;
+		dgsum[1] -= DOT(uv[0], vd) * gfact;
 	}
+	dg[0] = dgsum[0] / (hp->ns*hp->ns);
+	dg[1] = dgsum[1] / (hp->ns*hp->ns);
 }
 
 
@@ -524,29 +528,36 @@ doambient(				/* compute ambient component */
 		free(hp);
 		return(0);		/* no valid samples */
 	}
-	d = 1.0 / cnt;			/* final indirect irradiance/PI */
-	acol[0] *= d; acol[1] *= d; acol[2] *= d;
-	copycolor(rcol, acol);
+	copycolor(rcol, acol);		/* final indirect irradiance/PI */
 	if (cnt < hp->ns*hp->ns ||	/* incomplete sampling? */
 			(ra == NULL) & (pg == NULL) & (dg == NULL)) {
 		free(hp);
 		return(-1);		/* no radius or gradient calc. */
 	}
-	d = 0.01 * bright(rcol);	/* add in 1% before Hessian comp. */
-	if (d < FTINY) d = FTINY;
-	ap = hp->sa;			/* using Y channel from here on... */
+	multcolor(acol, hp->acoef);	/* normalize Y values */
+	if ((d = bright(acol)) > FTINY)
+		d = 1.0/d;
+	else
+		d = 0.0;
+	ap = hp->sa;			/* relative Y channel from here on... */
 	for (i = hp->ns*hp->ns; i--; ap++)
-		colval(ap->v,CIEY) = bright(ap->v) + d;
+		colval(ap->v,CIEY) = bright(ap->v)*d + 0.0314;
 
 	if (uv == NULL)			/* make sure we have axis pointers */
 		uv = my_uv;
 					/* compute radii & pos. gradient */
 	ambHessian(hp, uv, ra, pg);
+
 	if (dg != NULL)			/* compute direction gradient */
 		ambdirgrad(hp, uv, dg);
+
 	if (ra != NULL) {		/* scale/clamp radii */
-		d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt));
-		ra[0] *= d;
+		if (ra[0] < minarad) {
+			ra[0] = minarad;
+			if (ra[1] < minarad)
+				ra[1] = minarad;
+		}
+		ra[0] *= d = 1.0/sqrt(sqrt(wt));
 		if ((ra[1] *= d) > 2.0*ra[0])
 			ra[1] = 2.0*ra[0];
 		if (ra[1] > maxarad) {