--- ray/src/common/bsdf.c	2010/07/03 02:50:22	2.3
+++ ray/src/common/bsdf.c	2011/01/06 04:39:08	2.11
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf.c,v 2.3 2010/07/03 02:50:22 greg Exp $";
+static const char RCSid[] = "$Id: bsdf.c,v 2.11 2011/01/06 04:39:08 greg Exp $";
 #endif
 /*
  * Routines for handling BSDF data
@@ -23,7 +23,7 @@ typedef struct {
 	}	lat[MAXLATS+1];		/* latitudes */
 } ANGLE_BASIS;
 
-#define	MAXABASES	5		/* limit on defined bases */
+#define	MAXABASES	7		/* limit on defined bases */
 
 static ANGLE_BASIS	abase_list[MAXABASES] = {
 	{
@@ -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
@@ -200,7 +208,7 @@ ab_getndxR(		/* get index corresponding to the reverse
 
 
 static void
-load_angle_basis(	/* load BSDF angle basis */
+load_angle_basis(	/* load custom BSDF angle basis */
 	ezxml_t wab
 )
 {
@@ -208,11 +216,11 @@ load_angle_basis(	/* load BSDF angle basis */
 	ezxml_t	wbb;
 	int	i;
 	
-	if (abname == NULL || !*abname)
+	if (!abname || !*abname)
 		return;
 	for (i = nabases; i--; )
-		if (!strcmp(abname, abase_list[i].name))
-			return;		/* XXX assume it's the same */
+		if (!strcasecamp(abname, abase_list[i].name))
+			return;		/* assume it's the same */
 	if (nabases >= MAXABASES)
 		error(INTERNAL, "too many angle bases");
 	strcpy(abase_list[nabases].name, abname);
@@ -227,7 +235,7 @@ load_angle_basis(	/* load 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 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
@@ -250,13 +315,12 @@ load_bsdf_data(		/* load BSDF distribution for this wa
 	char  *sdata;
 	int  i;
 	
-	if ((cbasis == NULL || !*cbasis) | (rbasis == NULL || !*rbasis)) {
+	if ((!cbasis || !*cbasis) | (!rbasis || !*rbasis)) {
 		error(WARNING, "missing column/row basis for BSDF");
 		return;
 	}
-	/* XXX need to add routines for loading in foreign bases */
 	for (i = nabases; i--; )
-		if (!strcmp(cbasis, abase_list[i].name)) {
+		if (!strcasecamp(cbasis, abase_list[i].name)) {
 			dp->ninc = abase_list[i].nangles;
 			dp->ib_priv = (void *)&abase_list[i];
 			dp->ib_vec = ab_getvecR;
@@ -265,12 +329,12 @@ load_bsdf_data(		/* load BSDF distribution for this wa
 			break;
 		}
 	if (i < 0) {
-		sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
+		sprintf(errmsg, "undefined ColumnAngleBasis '%s'", cbasis);
 		error(WARNING, errmsg);
 		return;
 	}
 	for (i = nabases; i--; )
-		if (!strcmp(rbasis, abase_list[i].name)) {
+		if (!strcasecamp(rbasis, abase_list[i].name)) {
 			dp->nout = abase_list[i].nangles;
 			dp->ob_priv = (void *)&abase_list[i];
 			dp->ob_vec = ab_getvec;
@@ -279,13 +343,13 @@ load_bsdf_data(		/* load BSDF distribution for this wa
 			break;
 		}
 	if (i < 0) {
-		sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
+		sprintf(errmsg, "undefined RowAngleBasis '%s'", cbasis);
 		error(WARNING, errmsg);
 		return;
 	}
 				/* read BSDF data */
 	sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
-	if (sdata == NULL || !*sdata) {
+	if (!sdata || !*sdata) {
 		error(WARNING, "missing BSDF ScatteringData");
 		return;
 	}
@@ -309,7 +373,7 @@ load_bsdf_data(		/* load BSDF distribution for this wa
 		sdata++;
 	if (*sdata) {
 		sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
-				strlen(sdata));
+				(int)strlen(sdata));
 		error(WARNING, errmsg);
 	}
 }
@@ -321,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;
@@ -390,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);
@@ -401,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
@@ -454,16 +526,27 @@ 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 (strcasecamp(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")),
+		if (strcasecamp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
 					"Transmission Front"))
 			continue;
 		load_bsdf_data(dp, wdb);	/* load front BTDF */
@@ -487,6 +570,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);
@@ -593,12 +678,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) {
@@ -629,5 +716,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);
 }