--- ray/src/common/bsdf.c	2011/04/08 23:23:28	2.18
+++ ray/src/common/bsdf.c	2011/09/17 22:09:33	2.36
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf.c,v 2.18 2011/04/08 23:23:28 greg Exp $";
+static const char RCSid[] = "$Id: bsdf.c,v 2.36 2011/09/17 22:09:33 greg Exp $";
 #endif
 /*
  *  bsdf.c
@@ -10,9 +10,12 @@ static const char RCSid[] = "$Id: bsdf.c,v 2.18 2011/0
  *
  */
 
+#define	_USE_MATH_DEFINES
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include <math.h>
+#include <ctype.h>
 #include "ezxml.h"
 #include "hilbert.h"
 #include "bsdf.h"
@@ -45,10 +48,14 @@ int			SDretainSet = SDretainNone;
 SDError
 SDreportEnglish(SDError ec, FILE *fp)
 {
-	if (fp == NULL)
-		return ec;
 	if (!ec)
 		return SDEnone;
+	if ((ec < SDEnone) | (ec > SDEunknown)) {
+		SDerrorDetail[0] = '\0';
+		ec = SDEunknown;
+	}
+	if (fp == NULL)
+		return ec;
 	fputs(SDerrorEnglish[ec], fp);
 	if (SDerrorDetail[0]) {
 		fputs(": ", fp);
@@ -95,13 +102,17 @@ SDloadGeometry(SDData *sd, ezxml_t wdb)
 	if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
 		sd->dim[2] = atof(ezxml_txt(geom)) *
 				to_meters(ezxml_attr(geom, "unit"));
-	if ((sd->dim[0] < .0) | (sd->dim[1] < .0) | (sd->dim[2] < .0)) {
+	if ((sd->dim[0] < 0) | (sd->dim[1] < 0) | (sd->dim[2] < 0)) {
 		sprintf(SDerrorDetail, "Negative size in \"%s\"", sd->name);
 		return SDEdata;
 	}
 	if ((geom = ezxml_child(wdb, "Geometry")) == NULL ||
 			(mgfstr = ezxml_txt(geom)) == NULL)
 		return SDEnone;
+	while (isspace(*mgfstr))
+		++mgfstr;
+	if (!*mgfstr)
+		return SDEnone;
 	if ((fmt = ezxml_attr(geom, "format")) != NULL &&
 			strcasecmp(fmt, "MGF")) {
 		sprintf(SDerrorDetail,
@@ -160,20 +171,16 @@ SDloadFile(SDData *sd, const char *fname)
 	}
 				/* load geometry if present */
 	lastErr = SDloadGeometry(sd, ezxml_child(wtl, "Material"));
-	if (lastErr)
+	if (lastErr) {
+		ezxml_free(fl);
 		return lastErr;
+	}
 				/* try loading variable resolution data */
 	lastErr = SDloadTre(sd, wtl);
 				/* check our result */
-	switch (lastErr) {
-	case SDEformat:
-	case SDEdata:
-	case SDEsupport:	/* possibly we just tried the wrong format */
+	if (lastErr == SDEsupport)	/* try matrix BSDF if not tree data */
 		lastErr = SDloadMtx(sd, wtl);
-		break;
-	default:		/* variable res. OK else serious error */
-		break;
-	}
+		
 				/* done with XML file */
 	ezxml_free(fl);
 	
@@ -245,7 +252,8 @@ SDfreeSpectralDF(SDSpectralDF *df)
 		return;
 	SDfreeCumulativeCache(df);
 	for (n = df->ncomp; n-- > 0; )
-		(*df->comp[n].func->freeSC)(df->comp[n].dist);
+		if (df->comp[n].dist != NULL)
+			(*df->comp[n].func->freeSC)(df->comp[n].dist);
 	free(df);
 }
 
@@ -269,10 +277,14 @@ SDclipName(char *res, const char *fname)
 
 /* Initialize an unused BSDF struct (simply clears to zeroes) */
 void
-SDclearBSDF(SDData *sd)
+SDclearBSDF(SDData *sd, const char *fname)
 {
-	if (sd != NULL)
-		memset(sd, 0, sizeof(SDData));
+	if (sd == NULL)
+		return;
+	memset(sd, 0, sizeof(SDData));
+	if (fname == NULL)
+		return;
+	SDclipName(sd->name, fname);
 }
 
 /* Free data associated with BSDF struct */
@@ -330,7 +342,7 @@ SDgetCache(const char *bname)
 	sdl->next = SDcacheList;
 	SDcacheList = sdl;
 
-	sdl->refcnt++;
+	sdl->refcnt = 1;
 	return &sdl->bsdf;
 }
 
@@ -374,7 +386,7 @@ SDfreeCache(const SDData *sd)
 	for (sdl = SDcacheList; sdl != NULL; sdl = (sdLast=sdl)->next)
 		if (&sdl->bsdf == sd)
 			break;
-	if (sdl == NULL || --sdl->refcnt)
+	if (sdl == NULL || (sdl->refcnt -= (sdl->refcnt > 0)))
 		return;			/* missing or still in use */
 					/* keep unreferenced data? */
 	if (SDisLoaded(sd) && SDretainSet) {
@@ -397,34 +409,35 @@ SDfreeCache(const SDData *sd)
 
 /* Sample an individual BSDF component */
 SDError
-SDsampComponent(SDValue *sv, FVECT outVec, const FVECT inVec,
-			double randX, SDComponent *sdc)
+SDsampComponent(SDValue *sv, FVECT ioVec, double randX, SDComponent *sdc)
 {
 	float		coef[SDmaxCh];
 	SDError		ec;
+	FVECT		inVec;
 	const SDCDst	*cd;
 	double		d;
 	int		n;
 					/* check arguments */
-	if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL))
+	if ((sv == NULL) | (ioVec == NULL) | (sdc == NULL))
 		return SDEargument;
 					/* get cumulative distribution */
+	VCOPY(inVec, ioVec);
 	cd = (*sdc->func->getCDist)(inVec, sdc);
 	if (cd == NULL)
 		return SDEmemory;
-	if (cd->cTotal <= 1e-7) {	/* anything to sample? */
+	if (cd->cTotal <= 1e-6) {	/* anything to sample? */
 		sv->spec = c_dfcolor;
 		sv->cieY = .0;
-		memset(outVec, 0, 3*sizeof(double));
+		memset(ioVec, 0, 3*sizeof(double));
 		return SDEnone;
 	}
 	sv->cieY = cd->cTotal;
 					/* compute sample direction */
-	ec = (*sdc->func->sampCDist)(outVec, randX, cd);
+	ec = (*sdc->func->sampCDist)(ioVec, randX, cd);
 	if (ec)
 		return ec;
 					/* get BSDF color */
-	n = (*sdc->func->getBSDFs)(coef, outVec, inVec, sdc->dist);
+	n = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
 	if (n <= 0) {
 		strcpy(SDerrorDetail, "BSDF sample value error");
 		return SDEinternal;
@@ -451,7 +464,7 @@ SDmultiSamp(double t[], int n, double randX)
 	bitmask_t	ndx, coord[MS_MAXDIM];
 	
 	while (n > MS_MAXDIM)		/* punt for higher dimensions */
-		t[--n] = rand()*(1./RAND_MAX);
+		t[--n] = rand()*(1./(RAND_MAX+.5));
 	nBits = (8*sizeof(bitmask_t) - 1) / n;
 	ndx = randX * (double)((bitmask_t)1 << (nBits*n));
 					/* get coordinate on Hilbert curve */
@@ -459,7 +472,7 @@ SDmultiSamp(double t[], int n, double randX)
 					/* convert back to [0,1) range */
 	scale = 1. / (double)((bitmask_t)1 << nBits);
 	while (n--)
-		t[n] = scale * ((double)coord[n] + rand()*(1./RAND_MAX));
+		t[n] = scale * ((double)coord[n] + rand()*(1./(RAND_MAX+.5)));
 }
 
 #undef MS_MAXDIM
@@ -472,7 +485,7 @@ SDdiffuseSamp(FVECT outVec, int outFront, double randX
 	SDmultiSamp(outVec, 2, randX);
 	SDsquare2disk(outVec, outVec[0], outVec[1]);
 	outVec[2] = 1. - outVec[0]*outVec[0] - outVec[1]*outVec[1];
-	if (outVec[2] > .0)		/* a bit of paranoia */
+	if (outVec[2] > 0)		/* a bit of paranoia */
 		outVec[2] = sqrt(outVec[2]);
 	if (!outFront)			/* going out back? */
 		outVec[2] = -outVec[2];
@@ -480,13 +493,14 @@ SDdiffuseSamp(FVECT outVec, int outFront, double randX
 
 /* Query projected solid angle coverage for non-diffuse BSDF direction */
 SDError
-SDsizeBSDF(double *projSA, const FVECT vec, int qflags, const SDData *sd)
+SDsizeBSDF(double *projSA, const FVECT v1, const RREAL *v2,
+				int qflags, const SDData *sd)
 {
-	SDSpectralDF	*rdf;
+	SDSpectralDF	*rdf, *tdf;
 	SDError		ec;
 	int		i;
 					/* check arguments */
-	if ((projSA == NULL) | (vec == NULL) | (sd == NULL))
+	if ((projSA == NULL) | (v1 == NULL) | (sd == NULL))
 		return SDEargument;
 					/* initialize extrema */
 	switch (qflags) {
@@ -502,20 +516,26 @@ SDsizeBSDF(double *projSA, const FVECT vec, int qflags
 	case 0:
 		return SDEargument;
 	}
-	if (vec[2] > .0)		/* front surface query? */
+	if (v1[2] > 0)			/* front surface query? */
 		rdf = sd->rf;
 	else
 		rdf = sd->rb;
+	tdf = sd->tf;
+	if (v2 != NULL)			/* bidirectional? */
+		if (v1[2] > 0 ^ v2[2] > 0)
+			rdf = NULL;
+		else
+			tdf = NULL;
 	ec = SDEdata;			/* run through components */
 	for (i = (rdf==NULL) ? 0 : rdf->ncomp; i--; ) {
-		ec = (*rdf->comp[i].func->queryProjSA)(projSA, vec, qflags,
-							rdf->comp[i].dist);
+		ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2,
+						qflags, &rdf->comp[i]);
 		if (ec)
 			return ec;
 	}
-	for (i = (sd->tf==NULL) ? 0 : sd->tf->ncomp; i--; ) {
-		ec = (*sd->tf->comp[i].func->queryProjSA)(projSA, vec, qflags,
-							sd->tf->comp[i].dist);
+	for (i = (tdf==NULL) ? 0 : tdf->ncomp; i--; ) {
+		ec = (*tdf->comp[i].func->queryProjSA)(projSA, v1, v2,
+						qflags, &tdf->comp[i]);
 		if (ec)
 			return ec;
 	}
@@ -539,8 +559,8 @@ SDevalBSDF(SDValue *sv, const FVECT outVec, const FVEC
 	if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL))
 		return SDEargument;
 					/* whose side are we on? */
-	inFront = (inVec[2] > .0);
-	outFront = (outVec[2] > .0);
+	inFront = (inVec[2] > 0);
+	outFront = (outVec[2] > 0);
 					/* start with diffuse portion */
 	if (inFront & outFront) {
 		*sv = sd->rLambFront;
@@ -557,7 +577,7 @@ SDevalBSDF(SDValue *sv, const FVECT outVec, const FVEC
 	i = (sdf != NULL) ? sdf->ncomp : 0;
 	while (i-- > 0) {
 		nch = (*sdf->comp[i].func->getBSDFs)(coef, outVec, inVec,
-							sdf->comp[i].dist);
+							&sdf->comp[i]);
 		while (nch-- > 0) {
 			c_cmix(&sv->spec, sv->cieY, &sv->spec,
 					coef[nch], &sdf->comp[i].cspec[nch]);
@@ -581,7 +601,7 @@ SDdirectHemi(const FVECT inVec, int sflags, const SDDa
 	if ((inVec == NULL) | (sd == NULL))
 		return .0;
 					/* gather diffuse components */
-	if (inVec[2] > .0) {
+	if (inVec[2] > 0) {
 		hsum = sd->rLambFront.cieY;
 		rdf = sd->rf;
 	} else /* !inFront */ {
@@ -612,10 +632,10 @@ SDdirectHemi(const FVECT inVec, int sflags, const SDDa
 
 /* Sample BSDF direction based on the given random variable */
 SDError
-SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVec,
-			double randX, int sflags, const SDData *sd)
+SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int sflags, const SDData *sd)
 {
 	SDError		ec;
+	FVECT		inVec;
 	int		inFront;
 	SDSpectralDF	*rdf;
 	double		rdiff;
@@ -624,11 +644,12 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
 	SDComponent	*sdc;
 	const SDCDst	**cdarr = NULL;
 					/* check arguments */
-	if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL) |
-			(randX < .0) | (randX >= 1.))
+	if ((sv == NULL) | (ioVec == NULL) | (sd == NULL) |
+			(randX < 0) | (randX >= 1.))
 		return SDEargument;
 					/* whose side are we on? */
-	inFront = (inVec[2] > .0);
+	VCOPY(inVec, ioVec);
+	inFront = (inVec[2] > 0);
 					/* remember diffuse portions */
 	if (inFront) {
 		*sv = sd->rLambFront;
@@ -666,16 +687,16 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
 		}
 		sv->cieY += cdarr[i]->cTotal;
 	}
-	if (sv->cieY <= 1e-7) {		/* anything to sample? */
+	if (sv->cieY <= 1e-6) {		/* anything to sample? */
 		sv->cieY = .0;
-		memset(outVec, 0, 3*sizeof(double));
+		memset(ioVec, 0, 3*sizeof(double));
 		return SDEnone;
 	}
 					/* scale random variable */
 	randX *= sv->cieY;
 					/* diffuse reflection? */
 	if (randX < rdiff) {
-		SDdiffuseSamp(outVec, inFront, randX/rdiff);
+		SDdiffuseSamp(ioVec, inFront, randX/rdiff);
 		goto done;
 	}
 	randX -= rdiff;
@@ -683,7 +704,7 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
 	if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) {
 		if (randX < sd->tLamb.cieY) {
 			sv->spec = sd->tLamb.spec;
-			SDdiffuseSamp(outVec, !inFront, randX/sd->tLamb.cieY);
+			SDdiffuseSamp(ioVec, !inFront, randX/sd->tLamb.cieY);
 			goto done;
 		}
 		randX -= sd->tLamb.cieY;
@@ -695,11 +716,11 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
 		return SDEinternal;
 					/* compute sample direction */
 	sdc = (i < nr) ? &rdf->comp[i] : &sd->tf->comp[i-nr];
-	ec = (*sdc->func->sampCDist)(outVec, randX/cdarr[i]->cTotal, cdarr[i]);
+	ec = (*sdc->func->sampCDist)(ioVec, randX/cdarr[i]->cTotal, cdarr[i]);
 	if (ec)
 		return ec;
 					/* compute color */
-	j = (*sdc->func->getBSDFs)(coef, outVec, inVec, sdc->dist);
+	j = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
 	if (j <= 0) {
 		sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error",
 				sd->name);
@@ -726,10 +747,10 @@ SDcompXform(RREAL vMtx[3][3], const FVECT sNrm, const 
 	if ((vMtx == NULL) | (sNrm == NULL) | (uVec == NULL))
 		return SDEargument;
 	VCOPY(vMtx[2], sNrm);
-	if (normalize(vMtx[2]) == .0)
+	if (normalize(vMtx[2]) == 0)
 		return SDEargument;
 	fcross(vMtx[0], uVec, vMtx[2]);
-	if (normalize(vMtx[0]) == .0)
+	if (normalize(vMtx[0]) == 0)
 		return SDEargument;
 	fcross(vMtx[1], vMtx[2], vMtx[0]);
 	return SDEnone;
@@ -749,7 +770,7 @@ SDinvXform(RREAL iMtx[3][3], RREAL vMtx[3][3])
 	mTmp[0][1] = vMtx[2][1]*vMtx[0][2] - vMtx[2][2]*vMtx[0][1];
 	mTmp[0][2] = vMtx[1][2]*vMtx[0][1] - vMtx[1][1]*vMtx[0][2];
 	d = vMtx[0][0]*mTmp[0][0] + vMtx[1][0]*mTmp[0][1] + vMtx[2][0]*mTmp[0][2];
-	if (d == .0) {
+	if (d == 0) {
 		strcpy(SDerrorDetail, "Zero determinant in matrix inversion");
 		return SDEargument;
 	}
@@ -776,12 +797,12 @@ SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT i
 	if (vMtx == NULL) {		/* assume they just want to normalize */
 		if (resVec != inpVec)
 			VCOPY(resVec, inpVec);
-		return (normalize(resVec) > .0) ? SDEnone : SDEargument;
+		return (normalize(resVec) > 0) ? SDEnone : SDEargument;
 	}
 	vTmp[0] = DOT(vMtx[0], inpVec);
 	vTmp[1] = DOT(vMtx[1], inpVec);
 	vTmp[2] = DOT(vMtx[2], inpVec);
-	if (normalize(vTmp) == .0)
+	if (normalize(vTmp) == 0)
 		return SDEargument;
 	VCOPY(resVec, vTmp);
 	return SDEnone;
@@ -796,7 +817,6 @@ SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT i
 
 #include "standard.h"
 #include "paths.h"
-#include <ctype.h>
 
 #define MAXLATS		46		/* maximum number of latitudes */
 
@@ -853,7 +873,7 @@ static int	nabases = 3;	/* current number of defined b
 static int
 fequal(double a, double b)
 {
-	if (b != .0)
+	if (b != 0)
 		a = a/b - 1.;
 	return((a <= 1e-6) & (a >= -1e-6));
 }
@@ -916,7 +936,7 @@ ab_getndx(		/* get index corresponding to the given ve
 {
 	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
 	int	li, ndx;
-	double	pol, azi, d;
+	double	pol, azi;
 
 	if ((v[2] < -1.0) | (v[2] > 1.0))
 		return(-1);
@@ -1156,7 +1176,7 @@ check_bsdf_data(	/* check that BSDF data is sane */
 )
 {
 	double		*omega_iarr, *omega_oarr;
-	double		dom, contrib, hemi_total, full_total;
+	double		dom, hemi_total, full_total;
 	int		nneg;
 	FVECT		v;
 	int		i, o;
@@ -1171,7 +1191,7 @@ check_bsdf_data(	/* check that BSDF data is sane */
 	hemi_total = .0;
 	for (i = dp->ninc; i--; ) {
 		dom = getBSDF_incohm(dp,i);
-		if (dom <= .0) {
+		if (dom <= 0) {
 			error(WARNING, "zero/negative incoming solid angle");
 			continue;
 		}
@@ -1194,7 +1214,7 @@ check_bsdf_data(	/* check that BSDF data is sane */
 	hemi_total = .0;
 	for (o = dp->nout; o--; ) {
 		dom = getBSDF_outohm(dp,o);
-		if (dom <= .0) {
+		if (dom <= 0) {
 			error(WARNING, "zero/negative outgoing solid angle");
 			continue;
 		}
@@ -1218,7 +1238,7 @@ check_bsdf_data(	/* check that BSDF data is sane */
 		hemi_total = .0;
 		for (o = dp->nout; o--; ) {
 			double	f = BSDF_value(dp,i,o);
-			if (f >= .0)
+			if (f >= 0)
 				hemi_total += f*omega_oarr[o];
 			else {
 				nneg += (f < -FTINY);
@@ -1306,9 +1326,11 @@ load_BSDF(		/* load BSDF data from file */
 		error(WARNING, errmsg);
 		ezxml_free(fl);
 		return(NULL);
-	}		
-	load_angle_basis(ezxml_child(ezxml_child(wtl,
-				"DataDefinition"), "AngleBasis"));
+	}
+	for (wld = ezxml_child(ezxml_child(wtl,
+				"DataDefinition"), "AngleBasis");
+			wld != NULL; wld = wld->next)
+		load_angle_basis(wld);
 	dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
 	load_geometry(dp, ezxml_child(wtl, "Material"));
 	for (wld = ezxml_child(wtl, "WavelengthData");