--- ray/src/gen/mkillum2.c	2007/12/13 07:03:37	2.28
+++ ray/src/gen/mkillum2.c	2009/09/08 23:05:47	2.33
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: mkillum2.c,v 2.28 2007/12/13 07:03:37 greg Exp $";
+static const char	RCSid[] = "$Id: mkillum2.c,v 2.33 2009/09/08 23:05:47 greg Exp $";
 #endif
 /*
  * Routines to do the actual calculation for mkillum
@@ -12,9 +12,13 @@ static const char	RCSid[] = "$Id: mkillum2.c,v 2.28 20
 #include  "cone.h"
 #include  "source.h"
 
+#ifndef NBSDFSAMPS
+#define NBSDFSAMPS	256		/* BSDF resampling count */
+#endif
 
 COLORV *	distarr = NULL;		/* distribution array */
 int		distsiz = 0;
+COLORV *	direct_discount = NULL;	/* amount to take off direct */
 
 
 void
@@ -41,6 +45,27 @@ newdist(			/* allocate & clear distribution array */
 }
 
 
+static void
+new_discount()			/* allocate space for direct contrib. record */
+{
+	if (distsiz <= 0)
+		return;
+	direct_discount = (COLORV *)calloc(distsiz, sizeof(COLOR));
+	if (direct_discount == NULL)
+		error(SYSTEM, "out of memory in new_discount");
+}
+
+
+static void
+done_discount()			/* clear off direct contrib. record */
+{
+	if (direct_discount == NULL)
+		return;
+	free((void *)direct_discount);
+	direct_discount = NULL;
+}
+
+
 int
 process_ray(			/* process a ray result or report error */
 	RAY *r,
@@ -58,6 +83,11 @@ process_ray(			/* process a ray result or report error
 	multcolor(r->rcol, r->rcoef);	/* in case it's a source ray */
 	colp = &distarr[r->rno * 3];
 	addcolor(colp, r->rcol);
+	if (r->rsrc >= 0 &&		/* remember source contrib. */
+			direct_discount != NULL) {
+		colp = &direct_discount[r->rno * 3];
+		addcolor(colp, r->rcol);
+	}
 	return(1);
 }
 
@@ -127,6 +157,7 @@ srcsamps(			/* sample sources from this surface positi
 			d = 5.*FTINY;
 		for (i = 3; i--; )
 			sr.rorg[i] += d*nrm[i];
+		samplendx++;			/* increment sample counter */
 		if (!srcray(&sr, NULL, &si))
 			break;			/* end of sources */
 						/* index direction */
@@ -143,14 +174,9 @@ srcsamps(			/* sample sources from this surface positi
 		} else {
 			if (v[2] >= -FTINY)
 				continue;	/* only sample transmission */
-			d = 1.0 - v[2]*v[2];
-			i = d*nalt;
-			d = atan2(-v[1], -v[0])/(2.*PI);
-			if (d < 0.0) d += 1.0;
-			j = d*nazi + 0.5;
-			if (j >= nazi) j = 0;
-			sr.rno = i*nazi + j;
-			d = nalt*nazi/PI * -v[2];
+			v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2];
+			sr.rno = flatindex(v, nalt, nazi);
+			d = nalt*nazi*(1./PI) * v[2];
 		}
 		d *= si.dom;			/* solid angle correction */
 		scalecolor(sr.rcoef, d);
@@ -222,7 +248,26 @@ flatdir(		/* compute uniform hemispherical direction *
 	dv[2] = sqrt(1. - alt);
 }
 
+int
+flatindex(		/* compute index for hemispherical direction */
+	FVECT	dv,
+	int	nalt,
+	int	nazi
+)
+{
+	double	d;
+	int	i, j;
+	
+	d = 1.0 - dv[2]*dv[2];
+	i = d*nalt;
+	d = atan2(dv[1], dv[0]) * (0.5/PI);
+	if (d < 0.0) d += 1.0;
+	j = d*nazi + 0.5;
+	if (j >= nazi) j = 0;
+	return(i*nazi + j);
+}
 
+
 int
 my_default(	/* default illum action */
 	OBJREC  *ob,
@@ -350,7 +395,7 @@ my_face(		/* make an illum face */
 		    raysamp(dim[1], org, dir);
 		}
 				/* add in direct component? */
-	if (!directvis && il->flags & IL_LIGHT) {
+	if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
 		MAT4	ixfm;
 		if (il->sd == NULL) {
 			for (i = 3; i--; ) {
@@ -361,8 +406,13 @@ my_face(		/* make an illum face */
 			}
 			ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
 			ixfm[3][3] = 1.;
-		} else if (!invmat4(ixfm, xfm))
-			objerror(ob, INTERNAL, "cannot invert BSDF transform");
+		} else {
+			if (!invmat4(ixfm, xfm))
+				objerror(ob, INTERNAL,
+					"cannot invert BSDF transform");
+			if (!(il->flags & IL_LIGHT))
+				new_discount();
+		}
 		dim[0] = random();
 		nallow = 10*il->nsamps;
 		for (i = 0; i < il->nsamps; i++) {
@@ -392,6 +442,7 @@ my_face(		/* make an illum face */
 		nalt = sqrt(il->sd->nout/PI) + .5;
 		nazi = PI*nalt + .5;
 		redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm);
+		done_discount();
 		if (!il->sampdens)
 			il->sampdens = nalt*nazi/PI + .999;
 	}
@@ -552,7 +603,7 @@ my_ring(		/* make an illum ring */
 		    raysamp(dim[1], org, dir);
 		}
 				/* add in direct component? */
-	if (!directvis && il->flags & IL_LIGHT) {
+	if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
 		MAT4	ixfm;
 		if (il->sd == NULL) {
 			for (i = 3; i--; ) {
@@ -563,8 +614,13 @@ my_ring(		/* make an illum ring */
 			}
 			ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
 			ixfm[3][3] = 1.;
-		} else if (!invmat4(ixfm, xfm))
-			objerror(ob, INTERNAL, "cannot invert BSDF transform");
+		} else {
+			if (!invmat4(ixfm, xfm))
+				objerror(ob, INTERNAL,
+					"cannot invert BSDF transform");
+			if (!(il->flags & IL_LIGHT))
+				new_discount();
+		}
 		dim[0] = random();
 		for (i = 0; i < il->nsamps; i++) {
 					/* randomize location */
@@ -587,6 +643,7 @@ my_ring(		/* make an illum ring */
 		nalt = sqrt(il->sd->nout/PI) + .5;
 		nazi = PI*nalt + .5;
 		redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
+		done_discount();
 		if (!il->sampdens)
 			il->sampdens = nalt*nazi/PI + .999;
 	}
@@ -600,4 +657,110 @@ my_ring(		/* make an illum ring */
 				/* clean up */
 	freecone(ob);
 	return(1);
+}
+
+
+void
+redistribute(		/* pass distarr ray sums through BSDF */
+	struct BSDF_data *b,
+	int nalt,
+	int nazi,
+	FVECT u,
+	FVECT v,
+	FVECT w,
+	MAT4 xm
+)
+{
+	int	nout = 0;
+	MAT4	mymat, inmat;
+	COLORV	*idist;
+	COLORV	*cp;
+	FVECT	dv;
+	double	wt;
+	int	i, j, k, c, o;
+	COLOR	col, cinc;
+					/* copy incoming distribution */
+	if (b->ninc > distsiz)
+		error(INTERNAL, "error 1 in redistribute");
+	idist = (COLORV *)malloc(sizeof(COLOR)*b->ninc);
+	if (idist == NULL)
+		error(SYSTEM, "out of memory in redistribute");
+	memcpy(idist, distarr, sizeof(COLOR)*b->ninc);
+					/* compose direction transform */
+	for (i = 3; i--; ) {
+		mymat[i][0] = u[i];
+		mymat[i][1] = v[i];
+		mymat[i][2] = w[i];
+		mymat[i][3] = 0.;
+	}
+	mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
+	mymat[3][3] = 1.;
+	if (xm != NULL)
+		multmat4(mymat, xm, mymat);
+	for (i = 3; i--; ) {		/* make sure it's normalized */
+		wt = 1./sqrt(	mymat[0][i]*mymat[0][i] +
+				mymat[1][i]*mymat[1][i] +
+				mymat[2][i]*mymat[2][i]	);
+		for (j = 3; j--; )
+			mymat[j][i] *= wt;
+	}
+	if (!invmat4(inmat, mymat))	/* need inverse as well */
+		error(INTERNAL, "cannot invert BSDF transform");
+	newdist(nalt*nazi);		/* resample distribution */
+	for (i = b->ninc; i--; ) {
+		int	direct_out = -1;
+		COLOR	cdir;
+		getBSDF_incvec(dv, b, i);	/* compute incident irrad. */
+		multv3(dv, dv, mymat);
+		if (dv[2] < 0.0) {
+			dv[0] = -dv[0]; dv[1] = -dv[1]; dv[2] = -dv[2];
+			direct_out += (direct_discount != NULL);
+		}
+		wt = getBSDF_incohm(b, i);
+		wt *= dv[2];			/* solid_angle*cosine(theta) */
+		cp = &idist[3*i];
+		copycolor(cinc, cp);
+		scalecolor(cinc, wt);
+		if (!direct_out) {		/* discount direct contr. */
+			cp = &direct_discount[3*i];
+			copycolor(cdir, cp);
+			scalecolor(cdir, -wt);
+			if (b->nout != b->ninc)
+				direct_out = flatindex(dv, nalt, nazi);
+			else
+				direct_out = i;	/* assumes dist. mirroring */
+		}
+		for (k = nalt; k--; )		/* loop over distribution */
+		  for (j = nazi; j--; ) {
+		    int	rstart = random();
+		    for (c = NBSDFSAMPS; c--; ) {
+			double  sp[2];
+			multisamp(sp, 2, urand(rstart+c));
+			flatdir(dv, (k + sp[0])/nalt,
+					(j + .5 - sp[1])/nazi);
+			multv3(dv, dv, inmat);
+						/* evaluate BSDF @ outgoing */
+			o = getBSDF_outndx(b, dv);
+			if (o < 0) {
+				nout++;
+				continue;
+			}
+			wt = BSDF_value(b, i, o) * (1./NBSDFSAMPS);
+			copycolor(col, cinc);
+			if (b->nout != b->ninc)
+				o = k*nazi + j;
+			if (o == direct_out)
+				addcolor(col, cdir);	/* minus direct */
+			scalecolor(col, wt);
+			cp = &distarr[3*o];
+			addcolor(cp, col);	/* sum into distribution */
+		    }
+		  }
+	}
+	free(idist);			/* free temp space */
+	if (nout) {
+		sprintf(errmsg, "missing %.1f%% of BSDF directions",
+				100.*nout/(b->ninc*nalt*nazi*NBSDFSAMPS));
+		error(WARNING, errmsg);
+	}
 }