--- ray/src/gen/mkillum2.c	2009/03/04 00:12:25	2.29
+++ ray/src/gen/mkillum2.c	2011/08/16 18:09:53	2.37
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: mkillum2.c,v 2.29 2009/03/04 00:12:25 greg Exp $";
+static const char	RCSid[] = "$Id: mkillum2.c,v 2.37 2011/08/16 18:09:53 greg Exp $";
 #endif
 /*
  * Routines to do the actual calculation for mkillum
@@ -11,12 +11,17 @@ static const char	RCSid[] = "$Id: mkillum2.c,v 2.29 20
 #include  "face.h"
 #include  "cone.h"
 #include  "source.h"
+#include  "paths.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
 newdist(			/* allocate & clear distribution array */
 	int siz
@@ -24,14 +29,14 @@ newdist(			/* allocate & clear distribution array */
 {
 	if (siz <= 0) {
 		if (distsiz > 0)
-			free((void *)distarr);
+			free(distarr);
 		distarr = NULL;
 		distsiz = 0;
 		return;
 	}
 	if (distsiz < siz) {
 		if (distsiz > 0)
-			free((void *)distarr);
+			free(distarr);
 		distarr = (COLORV *)malloc(sizeof(COLOR)*siz);
 		if (distarr == NULL)
 			error(SYSTEM, "out of memory in newdist");
@@ -57,7 +62,7 @@ done_discount()			/* clear off direct contrib. record 
 {
 	if (direct_discount == NULL)
 		return;
-	free((void *)direct_discount);
+	free(direct_discount);
 	direct_discount = NULL;
 }
 
@@ -103,7 +108,7 @@ raysamp(			/* queue a ray sample */
 	VCOPY(myRay.rorg, org);
 	VCOPY(myRay.rdir, dir);
 	myRay.rmax = .0;
-	rayorigin(&myRay, PRIMARY, NULL, NULL);
+	rayorigin(&myRay, PRIMARY|SPECULAR, NULL, NULL);
 	myRay.rno = ndx;
 					/* queue ray, check result */
 	process_ray(&myRay, ray_pqueue(&myRay));
@@ -118,16 +123,14 @@ srcsamps(			/* sample sources from this surface positi
 	MAT4 ixfm
 )
 {
-	int  nalt, nazi;
+	int  nalt=1, nazi=1;
 	SRCINDEX  si;
 	RAY  sr;
 	FVECT	v;
 	double  d;
 	int  i, j;
 						/* get sampling density */
-	if (il->sampdens <= 0) {
-		nalt = nazi = 1;
-	} else {
+	if (il->sd == NULL && il->sampdens > 0) {
 		i = PI * il->sampdens;
 		nalt = sqrt(i/PI) + .5;
 		nazi = PI*nalt + .5;
@@ -151,8 +154,8 @@ srcsamps(			/* sample sources from this surface positi
 			d = -1.0001*il->thick - 5.*FTINY;
 		else
 			d = 5.*FTINY;
-		for (i = 3; i--; )
-			sr.rorg[i] += d*nrm[i];
+		VSUM(sr.rorg, sr.rorg, nrm, d);
+		samplendx++;			/* increment sample counter */
 		if (!srcray(&sr, NULL, &si))
 			break;			/* end of sources */
 						/* index direction */
@@ -243,6 +246,7 @@ flatdir(		/* compute uniform hemispherical direction *
 	dv[2] = sqrt(1. - alt);
 }
 
+
 int
 flatindex(		/* compute index for hemispherical direction */
 	FVECT	dv,
@@ -264,6 +268,51 @@ flatindex(		/* compute index for hemispherical directi
 
 
 int
+printgeom(		/* print out detailed geometry for BSDF */
+	struct BSDF_data *sd,
+	char *xfrot,
+	FVECT ctr,
+	double s1,
+	double s2
+)
+{
+        static char	mgftemp[] = TEMPLATE;
+	char		cmdbuf[64];
+	FILE		*fp;
+	double		sca;
+
+	if (sd == NULL || sd->mgf == NULL)
+		return(0);
+	if (sd->dim[0] <= FTINY || sd->dim[1] <= FTINY)
+		return(0);
+	if ((s1 > s2) ^ (sd->dim[0] > sd->dim[1])) {
+		sca = s1; s1 = s2; s2 = sca;
+	}
+	s1 /= sd->dim[0];
+	s2 /= sd->dim[1];
+	sca = s1 > s2 ? s1 : s2;
+	strcpy(mgftemp, TEMPLATE);
+	if ((fp = fopen(mktemp(mgftemp), "w")) == NULL)
+		error(SYSTEM, "cannot create temporary file for MGF");
+					/* prepend our transform */
+	fprintf(fp, "xf%s -s %.5f -t %.5g %.5g %.5g\n",
+			xfrot, sca, ctr[0], ctr[1], ctr[2]);
+					/* output given MGF description */
+	fputs(sd->mgf, fp);
+	fputs("\nxf\n", fp);
+	if (fclose(fp) == EOF)
+		error(SYSTEM, "error writing MGF temporary file");
+					/* execute mgf2rad to convert MGF */
+	strcpy(cmdbuf, "mgf2rad ");
+	strcpy(cmdbuf+8, mgftemp);
+	fflush(stdout);
+	system(cmdbuf);
+	unlink(mgftemp);		/* clean up */
+	return(1);
+}
+
+
+int
 my_default(	/* default illum action */
 	OBJREC  *ob,
 	struct illum_args  *il,
@@ -293,6 +342,7 @@ my_face(		/* make an illum face */
 	FVECT  u, v;
 	double  ur[2], vr[2];
 	MAT4  xfm;
+	char  xfrot[64];
 	int  nallow;
 	FACE  *fa;
 	int  i, j;
@@ -304,7 +354,7 @@ my_face(		/* make an illum face */
 	}
 				/* set up sampling */
 	if (il->sd != NULL) {
-		if (!getBSDF_xfm(xfm, fa->norm, il->udir)) {
+		if (!getBSDF_xfm(xfm, fa->norm, il->udir, xfrot)) {
 			objerror(ob, WARNING, "illegal up direction");
 			freeface(ob);
 			return(my_default(ob, il, nm));
@@ -346,6 +396,15 @@ my_face(		/* make an illum face */
 		if (r2 < vr[0]) vr[0] = r2;
 		if (r2 > vr[1]) vr[1] = r2;
 	}
+				/* output detailed geometry? */
+	if (!(il->flags & IL_LIGHT) && il->sd != NULL && il->sd->mgf != NULL &&
+			il->thick <= FTINY) {
+		for (j = 3; j--; )
+			org[j] = .5*(ur[0]+ur[1])*u[j] +
+					.5*(vr[0]+vr[1])*v[j] +
+					fa->offset*fa->norm[j];
+		printgeom(il->sd, xfrot, org, ur[1]-ur[0], vr[1]-vr[0]);
+	}
 	dim[0] = random();
 				/* sample polygon */
 	nallow = 5*n*il->nsamps;
@@ -390,7 +449,7 @@ my_face(		/* make an illum face */
 		    raysamp(dim[1], org, dir);
 		}
 				/* add in direct component? */
-	if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
+	if (il->flags & IL_LIGHT || il->sd != NULL) {
 		MAT4	ixfm;
 		if (il->sd == NULL) {
 			for (i = 3; i--; ) {
@@ -544,7 +603,7 @@ my_ring(		/* make an illum ring */
 	co = getcone(ob, 0);
 				/* set up sampling */
 	if (il->sd != NULL) {
-		if (!getBSDF_xfm(xfm, co->ad, il->udir)) {
+		if (!getBSDF_xfm(xfm, co->ad, il->udir, NULL)) {
 			objerror(ob, WARNING, "illegal up direction");
 			freecone(ob);
 			return(my_default(ob, il, nm));
@@ -598,7 +657,7 @@ my_ring(		/* make an illum ring */
 		    raysamp(dim[1], org, dir);
 		}
 				/* add in direct component? */
-	if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
+	if (il->flags & IL_LIGHT || il->sd != NULL) {
 		MAT4	ixfm;
 		if (il->sd == NULL) {
 			for (i = 3; i--; ) {
@@ -652,4 +711,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*(k*nazi + j)];
+			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);
+	}
 }