--- ray/src/rt/m_bsdf.c	2011/02/20 17:43:43	2.6
+++ ray/src/rt/m_bsdf.c	2012/09/02 15:33:15	2.22
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: m_bsdf.c,v 2.6 2011/02/20 17:43:43 greg Exp $";
+static const char RCSid[] = "$Id: m_bsdf.c,v 2.22 2012/09/02 15:33:15 greg Exp $";
 #endif
 /*
  *  Shading for materials with BSDFs taken from XML data files
@@ -50,7 +50,7 @@ static const char RCSid[] = "$Id: m_bsdf.c,v 2.6 2011/
  *  Arguments for MAT_BSDF are:
  *	6+	thick	BSDFfile	ux uy uz	funcfile	transform
  *	0
- *	0|3|9	rdf	gdf	bdf
+ *	0|3|6|9	rdf	gdf	bdf
  *		rdb	gdb	bdb
  *		rdt	gdt	bdt
  */
@@ -67,7 +67,7 @@ typedef struct {
 	RAY	*pr;		/* intersected ray */
 	FVECT	pnorm;		/* perturbed surface normal */
 	FVECT	vray;		/* local outgoing (return) vector */
-	double	sr_vpsa;	/* sqrt of BSDF projected solid angle */
+	double	sr_vpsa[2];	/* sqrt of BSDF projected solid angle extrema */
 	RREAL	toloc[3][3];	/* world to local BSDF coords */
 	RREAL	fromloc[3][3];	/* local BSDF coords to world */
 	double	thick;		/* surface thickness */
@@ -82,10 +82,8 @@ typedef struct {
 
 /* Jitter ray sample according to projected solid angle and specjitter */
 static void
-bsdf_jitter(FVECT vres, BSDFDAT *ndp)
+bsdf_jitter(FVECT vres, BSDFDAT *ndp, double sr_psa)
 {
-	double	sr_psa = ndp->sr_vpsa;
-
 	VCOPY(vres, ndp->vray);
 	if (specjitter < 1.)
 		sr_psa *= specjitter;
@@ -96,6 +94,72 @@ bsdf_jitter(FVECT vres, BSDFDAT *ndp)
 	normalize(vres);
 }
 
+/* Evaluate BSDF for direct component, returning true if OK to proceed */
+static int
+direct_bsdf_OK(COLOR cval, FVECT ldir, double omega, BSDFDAT *ndp)
+{
+	int	nsamp, ok = 0;
+	FVECT	vsrc, vsmp, vjit;
+	double	tomega;
+	double	sf, tsr, sd[2];
+	COLOR	csmp;
+	SDValue	sv;
+	SDError	ec;
+	int	i;
+					/* transform source direction */
+	if (SDmapDir(vsrc, ndp->toloc, ldir) != SDEnone)
+		return(0);
+					/* assign number of samples */
+	ec = SDsizeBSDF(&tomega, ndp->vray, vsrc, SDqueryMin, ndp->sd);
+	if (ec)
+		goto baderror;
+					/* check indirect over-counting */
+	if (ndp->thick != 0 && ndp->pr->crtype & (SPECULAR|AMBIENT)
+				&& vsrc[2] > 0 ^ ndp->vray[2] > 0) {
+		double	dx = vsrc[0] + ndp->vray[0];
+		double	dy = vsrc[1] + ndp->vray[1];
+		if (dx*dx + dy*dy <= omega+tomega)
+			return(0);
+	}
+	sf = specjitter * ndp->pr->rweight;
+	if (25.*tomega <= omega)
+		nsamp = 100.*sf + .5;
+	else
+		nsamp = 4.*sf*omega/tomega + .5;
+	nsamp += !nsamp;
+	setcolor(cval, .0, .0, .0);	/* sample our source area */
+	sf = sqrt(omega);
+	tsr = sqrt(tomega);
+	for (i = nsamp; i--; ) {
+		VCOPY(vsmp, vsrc);	/* jitter query directions */
+		if (nsamp > 1) {
+			multisamp(sd, 2, (i + frandom())/(double)nsamp);
+			vsmp[0] += (sd[0] - .5)*sf;
+			vsmp[1] += (sd[1] - .5)*sf;
+			if (normalize(vsmp) == 0) {
+				--nsamp;
+				continue;
+			}
+		}
+		bsdf_jitter(vjit, ndp, tsr);
+					/* compute BSDF */
+		ec = SDevalBSDF(&sv, vjit, vsmp, ndp->sd);
+		if (ec)
+			goto baderror;
+		if (sv.cieY <= FTINY)	/* worth using? */
+			continue;
+		cvt_sdcolor(csmp, &sv);
+		addcolor(cval, csmp);	/* average it in */
+		++ok;
+	}
+	sf = 1./(double)nsamp;
+	scalecolor(cval, sf);
+	return(ok);
+baderror:
+	objerror(ndp->mp, USER, transSDError(ec));
+	return(0);			/* gratis return */
+}
+
 /* Compute source contribution for BSDF (reflected & transmitted) */
 static void
 dir_bsdf(
@@ -106,10 +170,6 @@ dir_bsdf(
 )
 {
 	BSDFDAT		*np = (BSDFDAT *)nnp;
-	SDError		ec;
-	SDValue		sv;
-	FVECT		vsrc;
-	FVECT		vjit;
 	double		ldot;
 	double		dtmp;
 	COLOR		ctmp;
@@ -120,7 +180,7 @@ dir_bsdf(
 	if ((-FTINY <= ldot) & (ldot <= FTINY))
 		return;
 
-	if (ldot > .0 && bright(np->rdiff) > FTINY) {
+	if (ldot > 0 && bright(np->rdiff) > FTINY) {
 		/*
 		 *  Compute added diffuse reflected component.
 		 */
@@ -129,7 +189,7 @@ dir_bsdf(
 		scalecolor(ctmp, dtmp);
 		addcolor(cval, ctmp);
 	}
-	if (ldot < .0 && bright(np->tdiff) > FTINY) {
+	if (ldot < 0 && bright(np->tdiff) > FTINY) {
 		/*
 		 *  Compute added diffuse transmission.
 		 */
@@ -141,21 +201,14 @@ dir_bsdf(
 	/*
 	 *  Compute scattering coefficient using BSDF.
 	 */
-	if (SDmapDir(vsrc, np->toloc, ldir) != SDEnone)
+	if (!direct_bsdf_OK(ctmp, ldir, omega, np))
 		return;
-	bsdf_jitter(vjit, np);
-	ec = SDevalBSDF(&sv, vjit, vsrc, np->sd);
-	if (ec)
-		objerror(np->mp, USER, transSDError(ec));
-
-	if (sv.cieY <= FTINY)		/* not worth using? */
-		return;
-	cvt_sdcolor(ctmp, &sv);
-	if (ldot > .0) {		/* pattern only diffuse reflection */
+	if (ldot > 0) {		/* pattern only diffuse reflection */
 		COLOR	ctmp1, ctmp2;
-		dtmp = (np->pr->rod > .0) ? np->sd->rLambFront.cieY
+		dtmp = (np->pr->rod > 0) ? np->sd->rLambFront.cieY
 					: np->sd->rLambBack.cieY;
-		dtmp /= PI * sv.cieY;	/* diffuse fraction */
+					/* diffuse fraction */
+		dtmp /= PI * bright(ctmp);
 		copycolor(ctmp2, np->pr->pcol);
 		scalecolor(ctmp2, dtmp);
 		setcolor(ctmp1, 1.-dtmp, 1.-dtmp, 1.-dtmp);
@@ -180,10 +233,6 @@ dir_brdf(
 )
 {
 	BSDFDAT		*np = (BSDFDAT *)nnp;
-	SDError		ec;
-	SDValue		sv;
-	FVECT		vsrc;
-	FVECT		vjit;
 	double		ldot;
 	double		dtmp;
 	COLOR		ctmp, ctmp1, ctmp2;
@@ -207,20 +256,12 @@ dir_brdf(
 	/*
 	 *  Compute reflection coefficient using BSDF.
 	 */
-	if (SDmapDir(vsrc, np->toloc, ldir) != SDEnone)
+	if (!direct_bsdf_OK(ctmp, ldir, omega, np))
 		return;
-	bsdf_jitter(vjit, np);
-	ec = SDevalBSDF(&sv, vjit, vsrc, np->sd);
-	if (ec)
-		objerror(np->mp, USER, transSDError(ec));
-
-	if (sv.cieY <= FTINY)		/* not worth using? */
-		return;
-	cvt_sdcolor(ctmp, &sv);
 					/* pattern only diffuse reflection */
-	dtmp = (np->pr->rod > .0) ? np->sd->rLambFront.cieY
+	dtmp = (np->pr->rod > 0) ? np->sd->rLambFront.cieY
 				: np->sd->rLambBack.cieY;
-	dtmp /= PI * sv.cieY;		/* diffuse fraction */
+	dtmp /= PI * bright(ctmp);	/* diffuse fraction */
 	copycolor(ctmp2, np->pr->pcol);
 	scalecolor(ctmp2, dtmp);
 	setcolor(ctmp1, 1.-dtmp, 1.-dtmp, 1.-dtmp);
@@ -241,10 +282,6 @@ dir_btdf(
 )
 {
 	BSDFDAT		*np = (BSDFDAT *)nnp;
-	SDError		ec;
-	SDValue		sv;
-	FVECT		vsrc;
-	FVECT		vjit;
 	double		ldot;
 	double		dtmp;
 	COLOR		ctmp;
@@ -268,16 +305,8 @@ dir_btdf(
 	/*
 	 *  Compute scattering coefficient using BSDF.
 	 */
-	if (SDmapDir(vsrc, np->toloc, ldir) != SDEnone)
+	if (!direct_bsdf_OK(ctmp, ldir, omega, np))
 		return;
-	bsdf_jitter(vjit, np);
-	ec = SDevalBSDF(&sv, vjit, vsrc, np->sd);
-	if (ec)
-		objerror(np->mp, USER, transSDError(ec));
-
-	if (sv.cieY <= FTINY)		/* not worth using? */
-		return;
-	cvt_sdcolor(ctmp, &sv);
 					/* full pattern on transmission */
 	multcolor(ctmp, np->pr->pcol);
 	dtmp = -ldot * omega;
@@ -290,56 +319,53 @@ static int
 sample_sdcomp(BSDFDAT *ndp, SDComponent *dcp, int usepat)
 {
 	int	nstarget = 1;
-	int	nsent = 0;
+	int	nsent;
 	SDError	ec;
 	SDValue bsv;
-	double	sthick;
-	FVECT	vjit, vsmp;
+	double	xrand;
+	FVECT	vsmp;
 	RAY	sr;
-	int	ntrials;
 						/* multiple samples? */
 	if (specjitter > 1.5) {
 		nstarget = specjitter*ndp->pr->rweight + .5;
-		if (nstarget < 1)
-			nstarget = 1;
+		nstarget += !nstarget;
 	}
-						/* run through our trials */
-	for (ntrials = 0; nsent < nstarget && ntrials < 9*nstarget; ntrials++) {
-		SDerrorDetail[0] = '\0';
-						/* sample direction & coef. */
-		bsdf_jitter(vjit, ndp);
-		ec = SDsampComponent(&bsv, vsmp, vjit, ntrials ? frandom()
-				: urand(ilhash(dimlist,ndims)+samplendx), dcp);
+						/* run through our samples */
+	for (nsent = 0; nsent < nstarget; nsent++) {
+		if (nstarget == 1) {		/* stratify random variable */
+			xrand = urand(ilhash(dimlist,ndims)+samplendx);
+			if (specjitter < 1.)
+				xrand = .5 + specjitter*(xrand-.5);
+		} else {
+			xrand = (nsent + frandom())/(double)nstarget;
+		}
+		SDerrorDetail[0] = '\0';	/* sample direction & coef. */
+		bsdf_jitter(vsmp, ndp, ndp->sr_vpsa[0]);
+		ec = SDsampComponent(&bsv, vsmp, xrand, dcp);
 		if (ec)
 			objerror(ndp->mp, USER, transSDError(ec));
-						/* zero component? */
-		if (bsv.cieY <= FTINY)
+		if (bsv.cieY <= FTINY)		/* zero component? */
 			break;
 						/* map vector to world */
 		if (SDmapDir(sr.rdir, ndp->fromloc, vsmp) != SDEnone)
 			break;
-						/* unintentional penetration? */
-		if (DOT(sr.rdir, ndp->pr->ron) > .0 ^ vsmp[2] > .0)
-			continue;
 						/* spawn a specular ray */
 		if (nstarget > 1)
 			bsv.cieY /= (double)nstarget;
-		cvt_sdcolor(sr.rcoef, &bsv);	/* use color */
-		if (usepat)			/* pattern on transmission */
+		cvt_sdcolor(sr.rcoef, &bsv);	/* use sample color */
+		if (usepat)			/* apply pattern? */
 			multcolor(sr.rcoef, ndp->pr->pcol);
 		if (rayorigin(&sr, SPECULAR, ndp->pr, sr.rcoef) < 0) {
-			if (maxdepth  > 0)
+			if (maxdepth > 0)
 				break;
-			++nsent;		/* Russian roulette victim */
-			continue;
+			continue;		/* Russian roulette victim */
 		}
 						/* need to offset origin? */
-		if (ndp->thick != .0 && ndp->pr->rod > .0 ^ vsmp[2] > .0)
+		if (ndp->thick != 0 && ndp->pr->rod > 0 ^ vsmp[2] > 0)
 			VSUM(sr.rorg, sr.rorg, ndp->pr->ron, -ndp->thick);
 		rayvalue(&sr);			/* send & evaluate sample */
 		multcolor(sr.rcol, sr.rcoef);
 		addcolor(ndp->pr->rcol, sr.rcol);
-		++nsent;
 	}
 	return(nsent);
 }
@@ -354,11 +380,14 @@ sample_sdf(BSDFDAT *ndp, int sflags)
 
 	if (sflags == SDsampSpT) {
 		unsc = ndp->tunsamp;
-		dfp = ndp->sd->tf;
+		if (ndp->pr->rod > 0)
+			dfp = (ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb;
+		else
+			dfp = (ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf;
 		cvt_sdcolor(unsc, &ndp->sd->tLamb);
 	} else /* sflags == SDsampSpR */ {
 		unsc = ndp->runsamp;
-		if (ndp->pr->rod > .0) {
+		if (ndp->pr->rod > 0) {
 			dfp = ndp->sd->rf;
 			cvt_sdcolor(unsc, &ndp->sd->rLambFront);
 		} else {
@@ -375,7 +404,7 @@ sample_sdf(BSDFDAT *ndp, int sflags)
 			FVECT	vjit;
 			double	d;
 			COLOR	ctmp;
-			bsdf_jitter(vjit, ndp);
+			bsdf_jitter(vjit, ndp, ndp->sr_vpsa[1]);
 			d = SDdirectHemi(vjit, sflags, ndp->sd);
 			if (sflags == SDsampSpT) {
 				copycolor(ctmp, ndp->pr->pcol);
@@ -412,7 +441,7 @@ m_bsdf(OBJREC *m, RAY *r)
 				(m->oargs.nfargs % 3))
 		objerror(m, USER, "bad # arguments");
 						/* record surface struck */
-	hitfront = (r->rod > .0);
+	hitfront = (r->rod > 0);
 						/* load cal file */
 	mf = getfunc(m, 5, 0x1d, 1);
 						/* get thickness */
@@ -421,13 +450,13 @@ m_bsdf(OBJREC *m, RAY *r)
 		nd.thick = .0;
 						/* check shadow */
 	if (r->crtype & SHADOW) {
-		if (nd.thick != .0)
+		if (nd.thick != 0)
 			raytrans(r);		/* pass-through */
 		return(1);			/* or shadow */
 	}
 						/* check other rays to pass */
 	if (nd.thick != 0 && (!(r->crtype & (SPECULAR|AMBIENT)) ||
-				nd.thick > .0 ^ hitfront)) {
+				nd.thick > 0 ^ hitfront)) {
 		raytrans(r);			/* hide our proxy */
 		return(1);
 	}
@@ -444,7 +473,7 @@ m_bsdf(OBJREC *m, RAY *r)
 	} else {
 		if (m->oargs.nfargs < 6) {	/* check invisible backside */
 			if (!backvis && (nd.sd->rb == NULL) &
-						(nd.sd->tf == NULL)) {
+					(nd.sd->tb == NULL)) {
 				SDfreeCache(nd.sd);
 				raytrans(r);
 				return(1);
@@ -474,9 +503,9 @@ m_bsdf(OBJREC *m, RAY *r)
 	upvec[1] = evalue(mf->ep[2]);
 	upvec[2] = evalue(mf->ep[3]);
 						/* return to world coords */
-	if (mf->f != &unitxf) {
-		multv3(upvec, upvec, mf->f->xfm);
-		nd.thick *= mf->f->sca;
+	if (mf->fxp != &unitxf) {
+		multv3(upvec, upvec, mf->fxp->xfm);
+		nd.thick *= mf->fxp->sca;
 	}
 	raynormal(nd.pnorm, r);
 						/* compute local BSDF xform */
@@ -489,16 +518,17 @@ m_bsdf(OBJREC *m, RAY *r)
 	}
 	if (!ec)
 		ec = SDinvXform(nd.fromloc, nd.toloc);
-						/* determine BSDF resolution */
-	if (!ec)
-		ec = SDsizeBSDF(&nd.sr_vpsa, nd.vray, SDqueryMin, nd.sd);
-	if (!ec)
-		nd.sr_vpsa = sqrt(nd.sr_vpsa);
-	else {
-		objerror(m, WARNING, transSDError(ec));
-		SDfreeCache(nd.sd);
+	if (ec) {
+		objerror(m, WARNING, "Illegal orientation vector");
 		return(1);
 	}
+						/* determine BSDF resolution */
+	ec = SDsizeBSDF(nd.sr_vpsa, nd.vray, NULL, SDqueryMin+SDqueryMax, nd.sd);
+	if (ec)
+		objerror(m, USER, transSDError(ec));
+
+	nd.sr_vpsa[0] = sqrt(nd.sr_vpsa[0]);
+	nd.sr_vpsa[1] = sqrt(nd.sr_vpsa[1]);
 	if (!hitfront) {			/* perturb normal towards hit */
 		nd.pnorm[0] = -nd.pnorm[0];
 		nd.pnorm[1] = -nd.pnorm[1];
@@ -528,9 +558,9 @@ m_bsdf(OBJREC *m, RAY *r)
 		bnorm[0] = -nd.pnorm[0];
 		bnorm[1] = -nd.pnorm[1];
 		bnorm[2] = -nd.pnorm[2];
-		if (nd.thick != .0) {		/* proxy with offset? */
+		if (nd.thick != 0) {		/* proxy with offset? */
 			VCOPY(vtmp, r->rop);
-			VSUM(r->rop, vtmp, r->ron, -nd.thick);
+			VSUM(r->rop, vtmp, r->ron, nd.thick);
 			multambient(ctmp, r, bnorm);
 			VCOPY(r->rop, vtmp);
 		} else
@@ -540,9 +570,10 @@ m_bsdf(OBJREC *m, RAY *r)
 			flipsurface(r);
 	}
 						/* add direct component */
-	if ((bright(nd.tdiff) <= FTINY) & (nd.sd->tf == NULL)) {
+	if ((bright(nd.tdiff) <= FTINY) & (nd.sd->tf == NULL) &
+					(nd.sd->tb == NULL)) {
 		direct(r, dir_brdf, &nd);	/* reflection only */
-	} else if (nd.thick == .0) {
+	} else if (nd.thick == 0) {
 		direct(r, dir_bsdf, &nd);	/* thin surface scattering */
 	} else {
 		direct(r, dir_brdf, &nd);	/* reflection first */