--- ray/src/common/bsdf_m.c	2011/02/19 01:48:59	3.4
+++ ray/src/common/bsdf_m.c	2011/04/24 19:39:21	3.12
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf_m.c,v 3.4 2011/02/19 01:48:59 greg Exp $";
+static const char RCSid[] = "$Id: bsdf_m.c,v 3.12 2011/04/24 19:39:21 greg Exp $";
 #endif
 /*
  *  bsdf_m.c
@@ -81,12 +81,12 @@ 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);
 }
 
-/* returns the name of the given tag */
+/* Returns the name of the given tag */
 #ifdef ezxml_name
 #undef ezxml_name
 static char *
@@ -98,7 +98,7 @@ ezxml_name(ezxml_t xml)
 }
 #endif
 
-/* returns the given tag's character content or empty string if none */
+/* Returns the given tag's character content or empty string if none */
 #ifdef ezxml_txt
 #undef ezxml_txt
 static char *
@@ -157,16 +157,18 @@ SDnewMatrix(int ni, int no)
 /* Free a BSDF matrix */
 #define	SDfreeMatrix		free
 
-/* get vector for this angle basis index */
+/* get vector for this angle basis index (front exiting) */
 static int
-ab_getvec(FVECT v, int ndx, double randX, void *p)
+fo_getvec(FVECT v, double ndxr, void *p)
 {
-	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
+	ANGLE_BASIS	*ab = (ANGLE_BASIS *)p;
+	int		ndx = (int)ndxr;
+	double		randX = ndxr - ndx;
 	double		rx[2];
 	int		li;
 	double		pol, azi, d;
 	
-	if ((ndx < 0) | (ndx >= ab->nangles))
+	if ((ndxr < 0) | (ndx >= ab->nangles))
 		return RC_FAIL;
 	for (li = 0; ndx >= ab->lat[li].nphis; li++)
 		ndx -= ab->lat[li].nphis;
@@ -181,17 +183,17 @@ ab_getvec(FVECT v, int ndx, double randX, void *p)
 	return RC_GOOD;
 }
 
-/* get index corresponding to the given vector */
+/* get index corresponding to the given vector (front exiting) */
 static int
-ab_getndx(const FVECT v, void *p)
+fo_getndx(const FVECT v, void *p)
 {
-	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
+	ANGLE_BASIS	*ab = (ANGLE_BASIS *)p;
 	int	li, ndx;
 	double	pol, azi, d;
 
 	if (v == NULL)
 		return -1;
-	if ((v[2] < .0) | (v[2] > 1.0))
+	if ((v[2] < 0) | (v[2] > 1.))
 		return -1;
 	pol = 180.0/M_PI*acos(v[2]);
 	azi = 180.0/M_PI*atan2(v[1], v[0]);
@@ -210,9 +212,9 @@ ab_getndx(const FVECT v, void *p)
 /* compute square of real value */
 static double sq(double x) { return x*x; }
 
-/* get projected solid angle for this angle basis index */
+/* get projected solid angle for this angle basis index (universal) */
 static double
-ab_getohm(int ndx, void *p)
+io_getohm(int ndx, void *p)
 {
 	static int	last_li = -1;
 	static double	last_ohm;
@@ -235,13 +237,11 @@ ab_getohm(int ndx, void *p)
 				(double)ab->lat[li].nphis;
 }
 
-/* get reverse vector for this angle basis index */
+/* get vector for this angle basis index (back incident) */
 static int
-ab_getvecR(FVECT v, int ndx, double randX, void *p)
+bi_getvec(FVECT v, double ndxr, void *p)
 {
-	int	na = (*(ANGLE_BASIS *)p).nangles;
-
-	if (!ab_getvec(v, ndx, randX, p))
+	if (!fo_getvec(v, ndxr, p))
 		return RC_FAIL;
 
 	v[0] = -v[0];
@@ -251,9 +251,9 @@ ab_getvecR(FVECT v, int ndx, double randX, void *p)
 	return RC_GOOD;
 }
 
-/* get index corresponding to the reverse vector */
+/* get index corresponding to the vector (back incident) */
 static int
-ab_getndxR(const FVECT v, void *p)
+bi_getndx(const FVECT v, void *p)
 {
 	FVECT  v2;
 	
@@ -261,9 +261,60 @@ ab_getndxR(const FVECT v, void *p)
 	v2[1] = -v[1];
 	v2[2] = -v[2];
 
-	return ab_getndx(v2, p);
+	return fo_getndx(v2, p);
 }
 
+/* get vector for this angle basis index (back exiting) */
+static int
+bo_getvec(FVECT v, double ndxr, void *p)
+{
+	if (!fo_getvec(v, ndxr, p))
+		return RC_FAIL;
+
+	v[2] = -v[2];
+
+	return RC_GOOD;
+}
+
+/* get index corresponding to the vector (back exiting) */
+static int
+bo_getndx(const FVECT v, void *p)
+{
+	FVECT  v2;
+	
+	v2[0] = v[0];
+	v2[1] = v[1];
+	v2[2] = -v[2];
+
+	return fo_getndx(v2, p);
+}
+
+/* get vector for this angle basis index (front incident) */
+static int
+fi_getvec(FVECT v, double ndxr, void *p)
+{
+	if (!fo_getvec(v, ndxr, p))
+		return RC_FAIL;
+
+	v[0] = -v[0];
+	v[1] = -v[1];
+
+	return RC_GOOD;
+}
+
+/* get index corresponding to the vector (front incident) */
+static int
+fi_getndx(const FVECT v, void *p)
+{
+	FVECT  v2;
+	
+	v2[0] = -v[0];
+	v2[1] = -v[1];
+	v2[2] = v[2];
+
+	return fo_getndx(v2, p);
+}
+
 /* load custom BSDF angle basis */
 static int
 load_angle_basis(ezxml_t wab)
@@ -301,7 +352,7 @@ load_angle_basis(ezxml_t wab)
 					"ThetaBounds"), "LowerTheta"))),
 				abase_list[nabases].lat[i].tmin)) {
 			sprintf(SDerrorDetail, "Theta values disagree in '%s'",
-				abname);
+					abname);
 			return RC_DATERR;
 		}
 		abase_list[nabases].nangles +=
@@ -311,7 +362,7 @@ load_angle_basis(ezxml_t wab)
 				(abase_list[nabases].lat[i].nphis == 1 &&
 				abase_list[nabases].lat[i].tmin > FTINY)) {
 			sprintf(SDerrorDetail, "Illegal phi count in '%s'",
-				abname);
+					abname);
 			return RC_DATERR;
 		}
 	}
@@ -346,8 +397,7 @@ get_extrema(SDSpectralDF *df)
 	}
 	free(ohma);
 					/* need incoming solid angles, too? */
-	if (dp->ninc < dp->nout || dp->ib_ohm != dp->ob_ohm ||
-				dp->ib_priv != dp->ob_priv) {
+	if ((dp->ib_ohm != dp->ob_ohm) | (dp->ib_priv != dp->ob_priv)) {
 		double	ohm;
 		for (i = dp->ninc; i--; )
 			if ((ohm = mBSDF_incohm(dp,i)) < df->minProjSA)
@@ -363,28 +413,35 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
 	SDSpectralDF	*df;
 	SDMat		*dp;
 	char		*sdata;
+	int		tfront;
 	int		inbi, outbi;
 	int		i;
 					/* allocate BSDF component */
 	sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
-	if (!strcasecmp(sdata, "Transmission Front")) {
+	/*
+	 * Remember that front and back are reversed from WINDOW 6 orientations
+	 * Favor their "Front" (incoming light) since that's more often valid
+	 */
+	tfront = !strcasecmp(sdata, "Transmission Back");
+	if (!strcasecmp(sdata, "Transmission Front") ||
+			tfront & (sd->tf == NULL)) {
 		if (sd->tf != NULL)
 			SDfreeSpectralDF(sd->tf);
 		if ((sd->tf = SDnewSpectralDF(1)) == NULL)
 			return RC_MEMERR;
 		df = sd->tf;
 	} else if (!strcasecmp(sdata, "Reflection Front")) {
-		if (sd->rf != NULL)
-			SDfreeSpectralDF(sd->rf);
-		if ((sd->rf = SDnewSpectralDF(1)) == NULL)
-			return RC_MEMERR;
-		df = sd->rf;
-	} else if (!strcasecmp(sdata, "Reflection Back")) {
-		if (sd->rb != NULL)
+		if (sd->rb != NULL)	/* note back-front reversal */
 			SDfreeSpectralDF(sd->rb);
 		if ((sd->rb = SDnewSpectralDF(1)) == NULL)
 			return RC_MEMERR;
 		df = sd->rb;
+	} else if (!strcasecmp(sdata, "Reflection Back")) {
+		if (sd->rf != NULL)	/* note front-back reversal */
+			SDfreeSpectralDF(sd->rf);
+		if ((sd->rf = SDnewSpectralDF(1)) == NULL)
+			return RC_MEMERR;
+		df = sd->rf;
 	} else
 		return RC_FAIL;
 	/* XXX should also check "ScatteringDataType" for consistency? */
@@ -419,31 +476,38 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
 	dp = SDnewMatrix(abase_list[inbi].nangles, abase_list[outbi].nangles);
 	if (dp == NULL)
 		return RC_MEMERR;
-	dp->ib_priv = (void *)&abase_list[inbi];
-	dp->ob_priv = (void *)&abase_list[outbi];
+	dp->ib_priv = &abase_list[inbi];
+	dp->ob_priv = &abase_list[outbi];
 	if (df == sd->tf) {
-		dp->ib_vec = ab_getvecR;
-		dp->ib_ndx = ab_getndxR;
-		dp->ob_vec = ab_getvec;
-		dp->ob_ndx = ab_getndx;
+		if (tfront) {
+			dp->ib_vec = &fi_getvec;
+			dp->ib_ndx = &fi_getndx;
+			dp->ob_vec = &bo_getvec;
+			dp->ob_ndx = &bo_getndx;
+		} else {
+			dp->ib_vec = &bi_getvec;
+			dp->ib_ndx = &bi_getndx;
+			dp->ob_vec = &fo_getvec;
+			dp->ob_ndx = &fo_getndx;
+		}
 	} else if (df == sd->rf) {
-		dp->ib_vec = ab_getvec;
-		dp->ib_ndx = ab_getndx;
-		dp->ob_vec = ab_getvec;
-		dp->ob_ndx = ab_getndx;
+		dp->ib_vec = &fi_getvec;
+		dp->ib_ndx = &fi_getndx;
+		dp->ob_vec = &fo_getvec;
+		dp->ob_ndx = &fo_getndx;
 	} else /* df == sd->rb */ {
-		dp->ib_vec = ab_getvecR;
-		dp->ib_ndx = ab_getndxR;
-		dp->ob_vec = ab_getvecR;
-		dp->ob_ndx = ab_getndxR;
+		dp->ib_vec = &bi_getvec;
+		dp->ib_ndx = &bi_getndx;
+		dp->ob_vec = &bo_getvec;
+		dp->ob_ndx = &bo_getndx;
 	}
-	dp->ib_ohm = ab_getohm;
-	dp->ob_ohm = ab_getohm;
+	dp->ib_ohm = &io_getohm;
+	dp->ob_ohm = &io_getohm;
 	df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
 	df->comp[0].dist = dp;
 	df->comp[0].func = &SDhandleMtx;
 					/* read BSDF data */
-	sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
+	sdata = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
 	if (!sdata || !*sdata) {
 		sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
 				sd->name);
@@ -567,10 +631,14 @@ SDloadMtx(SDData *sd, ezxml_t wtl)
 /* Get Matrix BSDF value */
 static int
 SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
-				const FVECT inVec, const void *dist)
+				const FVECT inVec, SDComponent *sdc)
 {
-	const SDMat	*dp = (const SDMat *)dist;
+	const SDMat	*dp;
 	int		i_ndx, o_ndx;
+					/* check arguments */
+	if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
+			|| (dp = (SDMat *)sdc->dist) == NULL)
+		return 0;
 					/* get angle indices */
 	i_ndx = mBSDF_incndx(dp, inVec);
 	o_ndx = mBSDF_outndx(dp, outVec);
@@ -585,59 +653,52 @@ SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
 	return 1;			/* XXX monochrome for now */
 }
 
-/* Query solid angle for vector */
+/* Query solid angle for vector(s) */
 static SDError
-SDqueryMtxProjSA(double *psa, const FVECT vec, int qflags, const void *dist)
+SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2,
+					int qflags, SDComponent *sdc)
 {
-	const SDMat	*dp = (const SDMat *)dist;
-
-	if (!(qflags & SDqueryInc+SDqueryOut))
+	const SDMat	*dp;
+	double		inc_psa, out_psa;
+					/* check arguments */
+	if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
+			(dp = (SDMat *)sdc->dist) == NULL)
 		return SDEargument;
-	if (qflags & SDqueryInc) {
-		double	inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec));
-		if (inc_psa < .0)
-			return SDEinternal;
-		switch (qflags & SDqueryMin+SDqueryMax) {
-		case SDqueryMax:
-			if (inc_psa > psa[0])
-				psa[0] = inc_psa;
-			break;
-		case SDqueryMin+SDqueryMax:
-			if (inc_psa > psa[1])
-				psa[1] = inc_psa;
-			/* fall through */
-		case SDqueryMin:
-			if (inc_psa < psa[0])
-				psa[0] = inc_psa;
-			break;
-		case 0:
+	if (v2 == NULL)
+		v2 = v1;
+					/* get projected solid angles */
+	out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1));
+	inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2));
+	if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) {
+		inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2));
+		out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1));
+	}
+
+	switch (qflags) {		/* record based on flag settings */
+	case SDqueryMax:
+		if (inc_psa > psa[0])
 			psa[0] = inc_psa;
-			break;
-		}
+		if (out_psa > psa[0])
+			psa[0] = out_psa;
+		break;
+	case SDqueryMin+SDqueryMax:
+		if (inc_psa > psa[0])
+			psa[1] = inc_psa;
+		if (out_psa > psa[0])
+			psa[1] = out_psa;
+		/* fall through */
+	case SDqueryMin:
+	case SDqueryVal:
+		if (qflags == SDqueryVal)
+			psa[0] = M_PI;
+		if ((inc_psa > 0) & (inc_psa < psa[0]))
+			psa[0] = inc_psa;
+		if ((out_psa > 0) & (out_psa < psa[0]))
+			psa[0] = out_psa;
+		break;
 	}
-	if (qflags & SDqueryOut) {
-		double	out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec));
-		if (out_psa < .0)
-			return SDEinternal;
-		switch (qflags & SDqueryMin+SDqueryMax) {
-		case SDqueryMax:
-			if (out_psa > psa[0])
-				psa[0] = out_psa;
-			break;
-		case SDqueryMin+SDqueryMax:
-			if (out_psa > psa[1])
-				psa[1] = out_psa;
-			/* fall through */
-		case SDqueryMin:
-			if (out_psa < psa[0])
-				psa[0] = out_psa;
-			break;
-		case 0:
-			psa[(qflags&SDqueryInc)!=0] = out_psa;
-			break;
-		}
-	}
-	return SDEnone;
+					/* make sure it's legal */
+	return (psa[0] <= 0) ? SDEinternal : SDEnone;
 }
 
 /* Compute new cumulative distribution from BSDF */
@@ -675,12 +736,13 @@ make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp
 static const SDCDst *
 SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 {
-	SDMat		*dp = (SDMat *)sdc->dist;
+	SDMat		*dp;
 	int		reverse;
 	SDMatCDst	myCD;
 	SDMatCDst	*cd, *cdlast;
-
-	if (dp == NULL)
+					/* check arguments */
+	if ((inVec == NULL) | (sdc == NULL) ||
+			(dp = (SDMat *)sdc->dist) == NULL)
 		return NULL;
 	memset(&myCD, 0, sizeof(myCD));
 	myCD.indx = mBSDF_incndx(dp, inVec);
@@ -729,12 +791,15 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 
 /* Sample cumulative distribution */
 static SDError
-SDsampMtxCDist(FVECT outVec, double randX, const SDCDst *cdp)
+SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp)
 {
 	const unsigned	maxval = ~0;
 	const SDMatCDst	*mcd = (const SDMatCDst *)cdp;
 	const unsigned	target = randX*maxval;
 	int		i, iupper, ilower;
+					/* check arguments */
+	if ((ioVec == NULL) | (mcd == NULL))
+		return SDEargument;
 					/* binary search to find index */
 	ilower = 0; iupper = mcd->calen;
 	while ((i = (iupper + ilower) >> 1) != ilower)
@@ -746,9 +811,9 @@ SDsampMtxCDist(FVECT outVec, double randX, const SDCDs
 	randX = (randX*maxval - mcd->carr[ilower]) /
 			(double)(mcd->carr[iupper] - mcd->carr[ilower]);
 					/* convert index to vector */
-	if ((*mcd->ob_vec)(outVec, i, randX, mcd->ob_priv))
+	if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv))
 		return SDEnone;
-	strcpy(SDerrorDetail, "BSDF sampling fault");
+	strcpy(SDerrorDetail, "Matrix BSDF sampling fault");
 	return SDEinternal;
 }