--- ray/src/common/bsdf.c	2010/08/01 22:26:37	2.5
+++ ray/src/common/bsdf.c	2011/01/06 22:27:56	2.13
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf.c,v 2.5 2010/08/01 22:26:37 greg Exp $";
+static const char RCSid[] = "$Id: bsdf.c,v 2.13 2011/01/06 22:27:56 greg Exp $";
 #endif
 /*
  * Routines for handling BSDF data
@@ -61,8 +61,16 @@ static ANGLE_BASIS	abase_list[MAXABASES] = {
 
 static int	nabases = 3;	/* current number of defined bases */
 
-#define  FEQ(a,b)	((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
+#define  FEQ(a,b)	((a)-(b) <= 1e-6 && (b)-(a) <= 1e-6)
 
+static int
+fequal(double a, double b)
+{
+	if (b != .0)
+		a = a/b - 1.;
+	return((a <= 1e-6) & (a >= -1e-6));
+}
+
 // returns the name of the given tag
 #ifdef ezxml_name
 #undef ezxml_name
@@ -211,7 +219,7 @@ load_angle_basis(	/* load custom BSDF angle basis */
 	if (!abname || !*abname)
 		return;
 	for (i = nabases; i--; )
-		if (!strcmp(abname, abase_list[i].name))
+		if (!strcasecmp(abname, abase_list[i].name))
 			return;		/* assume it's the same */
 	if (nabases >= MAXABASES)
 		error(INTERNAL, "too many angle bases");
@@ -227,7 +235,7 @@ load_angle_basis(	/* load custom BSDF angle basis */
 		if (!i)
 			abase_list[nabases].lat[i].tmin =
 					-abase_list[nabases].lat[i+1].tmin;
-		else if (!FEQ(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
+		else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
 					"ThetaBounds"), "LowerTheta"))),
 				abase_list[nabases].lat[i].tmin))
 			error(WARNING, "theta values disagree in custom basis");
@@ -239,7 +247,64 @@ load_angle_basis(	/* load custom BSDF angle basis */
 }
 
 
+static double
+to_meters(		/* return factor to convert given unit to meters */
+	const char *unit
+)
+{
+	if (unit == NULL) return(1.);		/* safe assumption? */
+	if (!strcasecmp(unit, "Meter")) return(1.);
+	if (!strcasecmp(unit, "Foot")) return(.3048);
+	if (!strcasecmp(unit, "Inch")) return(.0254);
+	if (!strcasecmp(unit, "Centimeter")) return(.01);
+	if (!strcasecmp(unit, "Millimeter")) return(.001);
+	sprintf(errmsg, "unknown dimensional unit '%s'", unit);
+	error(USER, errmsg);
+}
+
+
 static void
+load_geometry(		/* load geometric dimensions and description (if any) */
+	struct BSDF_data *dp,
+	ezxml_t wdb
+)
+{
+	ezxml_t		geom;
+	double		cfact;
+	const char	*fmt, *mgfstr;
+
+	dp->dim[0] = dp->dim[1] = dp->dim[2] = 0;
+	dp->mgf = NULL;
+	if ((geom = ezxml_child(wdb, "Width")) != NULL)
+		dp->dim[0] = atof(ezxml_txt(geom)) *
+				to_meters(ezxml_attr(geom, "unit"));
+	if ((geom = ezxml_child(wdb, "Height")) != NULL)
+		dp->dim[1] = atof(ezxml_txt(geom)) *
+				to_meters(ezxml_attr(geom, "unit"));
+	if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
+		dp->dim[2] = atof(ezxml_txt(geom)) *
+				to_meters(ezxml_attr(geom, "unit"));
+	if ((geom = ezxml_child(wdb, "Geometry")) == NULL ||
+			(mgfstr = ezxml_txt(geom)) == NULL)
+		return;
+	if ((fmt = ezxml_attr(geom, "format")) != NULL &&
+			strcasecmp(fmt, "MGF")) {
+		sprintf(errmsg, "unrecognized geometry format '%s'", fmt);
+		error(WARNING, errmsg);
+		return;
+	}
+	cfact = to_meters(ezxml_attr(geom, "unit"));
+	dp->mgf = (char *)malloc(strlen(mgfstr)+32);
+	if (dp->mgf == NULL)
+		error(SYSTEM, "out of memory in load_geometry");
+	if (cfact < 0.99 || cfact > 1.01)
+		sprintf(dp->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
+	else
+		strcpy(dp->mgf, mgfstr);
+}
+
+
+static void
 load_bsdf_data(		/* load BSDF distribution for this wavelength */
 	struct BSDF_data *dp,
 	ezxml_t wdb
@@ -255,7 +320,7 @@ load_bsdf_data(		/* load BSDF distribution for this wa
 		return;
 	}
 	for (i = nabases; i--; )
-		if (!strcmp(cbasis, abase_list[i].name)) {
+		if (!strcasecmp(cbasis, abase_list[i].name)) {
 			dp->ninc = abase_list[i].nangles;
 			dp->ib_priv = (void *)&abase_list[i];
 			dp->ib_vec = ab_getvecR;
@@ -269,7 +334,7 @@ load_bsdf_data(		/* load BSDF distribution for this wa
 		return;
 	}
 	for (i = nabases; i--; )
-		if (!strcmp(rbasis, abase_list[i].name)) {
+		if (!strcasecmp(rbasis, abase_list[i].name)) {
 			dp->nout = abase_list[i].nangles;
 			dp->ob_priv = (void *)&abase_list[i];
 			dp->ob_vec = ab_getvec;
@@ -320,7 +385,7 @@ check_bsdf_data(	/* check that BSDF data is sane */
 )
 {
 	double		*omega_iarr, *omega_oarr;
-	double		dom, contrib, hemi_total;
+	double		dom, contrib, hemi_total, full_total;
 	int		nneg;
 	FVECT		v;
 	int		i, o;
@@ -389,7 +454,7 @@ check_bsdf_data(	/* check that BSDF data is sane */
 				BSDF_value(dp,i,o) = .0f;
 			}
 		}
-		if (hemi_total > 1.02) {
+		if (hemi_total > 1.01) {
 			sprintf(errmsg,
 			"incoming BSDF direction %d passes %.1f%% of light",
 					i, 100.*hemi_total);
@@ -400,23 +465,31 @@ check_bsdf_data(	/* check that BSDF data is sane */
 		sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
 		error(WARNING, errmsg);
 	}
-					/* reverse roles and check again */
+	full_total = .0;		/* reverse roles and check again */
 	for (o = 0; o < dp->nout; o++) {
 		hemi_total = .0;
 		for (i = dp->ninc; i--; )
 			hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
 
-		if (hemi_total > 1.02) {
+		if (hemi_total > 1.01) {
 			sprintf(errmsg,
 			"outgoing BSDF direction %d collects %.1f%% of light",
 					o, 100.*hemi_total);
 			error(WARNING, errmsg);
 		}
+		full_total += hemi_total*omega_oarr[o];
 	}
+	full_total /= PI;
+	if (full_total > 1.00001) {
+		sprintf(errmsg, "BSDF transfers %.4f%% of light",
+				100.*full_total);
+		error(WARNING, errmsg);
+	}
 	free(omega_iarr); free(omega_oarr);
 	return(1);
 }
 
+
 struct BSDF_data *
 load_BSDF(		/* load BSDF data from file */
 	char *fname
@@ -453,20 +526,35 @@ load_BSDF(		/* load BSDF data from file */
 		return(NULL);
 	}
 	wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
+	if (strcasecmp(ezxml_txt(ezxml_child(ezxml_child(wtl,
+			"DataDefinition"), "IncidentDataStructure")),
+			"Columns")) {
+		sprintf(errmsg,
+			"BSDF \"%s\": unsupported IncidentDataStructure",
+				path);
+		error(WARNING, errmsg);
+		ezxml_free(fl);
+		return(NULL);
+	}		
 	load_angle_basis(ezxml_child(ezxml_child(wtl,
 				"DataDefinition"), "AngleBasis"));
 	dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
+	load_geometry(dp, ezxml_child(wtl, "Material"));
 	for (wld = ezxml_child(wtl, "WavelengthData");
 				wld != NULL; wld = wld->next) {
-		if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
+		if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
+				"Visible"))
 			continue;
-		wdb = ezxml_child(wld, "WavelengthDataBlock");
-		if (wdb == NULL) continue;
-		if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
+		for (wdb = ezxml_child(wld, "WavelengthDataBlock");
+				wdb != NULL; wdb = wdb->next)
+			if (!strcasecmp(ezxml_txt(ezxml_child(wdb,
+					"WavelengthDataDirection")),
 					"Transmission Front"))
-			continue;
-		load_bsdf_data(dp, wdb);	/* load front BTDF */
-		break;				/* ignore the rest */
+				break;
+		if (wdb != NULL) {		/* load front BTDF */
+			load_bsdf_data(dp, wdb);
+			break;			/* ignore the rest */
+		}
 	}
 	ezxml_free(fl);				/* done with XML file */
 	if (!check_bsdf_data(dp)) {
@@ -486,6 +574,8 @@ free_BSDF(		/* free BSDF data structure */
 {
 	if (b == NULL)
 		return;
+	if (b->mgf != NULL)
+		free(b->mgf);
 	if (b->bsdf != NULL)
 		free(b->bsdf);
 	free(b);
@@ -592,12 +682,14 @@ int
 getBSDF_xfm(		/* compute BSDF orient. -> world orient. transform */
 	MAT4 xm,
 	FVECT nrm,
-	UpDir ud
+	UpDir ud,
+	char *xfbuf
 )
 {
 	char	*xfargs[7];
 	XF	myxf;
 	FVECT	updir, xdest, ydest;
+	int	i;
 
 	updir[0] = updir[1] = updir[2] = 0.;
 	switch (ud) {
@@ -628,5 +720,13 @@ getBSDF_xfm(		/* compute BSDF orient. -> world orient.
 	fcross(ydest, nrm, xdest);
 	xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
 	copymat4(xm, myxf.xfm);
+	if (xfbuf == NULL)
+		return(1);
+				/* return xf arguments as well */
+	for (i = 0; xfargs[i] != NULL; i++) {
+		*xfbuf++ = ' ';
+		strcpy(xfbuf, xfargs[i]);
+		while (*xfbuf) ++xfbuf;
+	}
 	return(1);
 }