--- ray/src/cv/bsdf2rad.c	2013/10/31 18:03:13	2.2
+++ ray/src/cv/bsdf2rad.c	2014/02/20 09:21:32	2.7
@@ -1,14 +1,16 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf2rad.c,v 2.2 2013/10/31 18:03:13 greg Exp $";
+static const char RCSid[] = "$Id: bsdf2rad.c,v 2.7 2014/02/20 09:21:32 greg Exp $";
 #endif
 /*
- *  Plot 3-D BSDF output based on scattering interpolant representation
+ *  Plot 3-D BSDF output based on scattering interpolant or XML representation
  */
 
 #define _USE_MATH_DEFINES
 #include <stdio.h>
+#include <string.h>
 #include <stdlib.h>
 #include <math.h>
+#include "rtprocess.h"
 #include "bsdfrep.h"
 
 const float	colarr[6][3] = {
@@ -27,56 +29,111 @@ int
 main(int argc, char *argv[])
 {
 	int	showPeaks = 0;
-	char	buf[128];
+	int	doTrans = 0;
+	int	inpXML = -1;
+	RBFNODE	*rbf = NULL;
 	FILE	*fp;
-	RBFNODE	*rbf;
+	char	buf[128];
+	SDData	myBSDF;
 	double	bsdf, min_log;
-	FVECT	dir;
+	FVECT	idir, odir;
 	int	i, j, n;
-
+						/* check arguments */
 	progname = argv[0];
 	if (argc > 1 && !strcmp(argv[1], "-p")) {
 		++showPeaks;
 		++argv; --argc;
 	}
-	if (argc < 4) {
+	if (argc > 1 && !strcmp(argv[1], "-t")) {
+		++doTrans;
+		++argv; --argc;
+	}
+	if (argc >= 4 && (n = strlen(argv[1])-4) > 0) {
+		if (!strcasecmp(argv[1]+n, ".xml"))
+			inpXML = 1;
+		else if (!strcasecmp(argv[1]+n, ".sir"))
+			inpXML = 0;
+	}
+	if (inpXML < 0) {
 		fprintf(stderr, "Usage: %s [-p] bsdf.sir theta1 phi1 .. > output.rad\n", progname);
+		fprintf(stderr, "   Or: %s [-t] bsdf.xml theta1 phi1 .. > output.rad\n", progname);
 		return(1);
 	}
 						/* load input */
-	if ((fp = fopen(argv[1], "rb")) == NULL) {
-		fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
-				progname, argv[1]);
-		return(1);
+	if (inpXML) {
+		SDclearBSDF(&myBSDF, argv[1]);
+		if (SDreportError(SDloadFile(&myBSDF, argv[1]), stderr))
+			return(1);
+		bsdf_min = 1./M_PI;
+		if (myBSDF.rf != NULL && myBSDF.rLambFront.cieY < bsdf_min*M_PI)
+			bsdf_min = myBSDF.rLambFront.cieY/M_PI;
+		if (myBSDF.rb != NULL && myBSDF.rLambBack.cieY < bsdf_min*M_PI)
+			bsdf_min = myBSDF.rLambBack.cieY/M_PI;
+		if ((myBSDF.tf != NULL) | (myBSDF.tb != NULL) &&
+				myBSDF.tLamb.cieY < bsdf_min*M_PI)
+			bsdf_min = myBSDF.tLamb.cieY/M_PI;
+		if (doTrans && (myBSDF.tf == NULL) & (myBSDF.tb == NULL)) {
+			fprintf(stderr, "%s: no transmitted component in '%s'\n",
+					progname, argv[1]);
+			return(1);
+		}
+	} else {
+		fp = fopen(argv[1], "rb");
+		if (fp == NULL) {
+			fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
+					progname, argv[1]);
+			return(1);
+		}
+		if (!load_bsdf_rep(fp))
+			return(1);
+		fclose(fp);
 	}
-	if (!load_bsdf_rep(fp))
-		return(1);
-	fclose(fp);
-	min_log = log(bsdf_min*.5);
+#ifdef DEBUG
+	fprintf(stderr, "Minimum BSDF set to %.4f\n", bsdf_min);
+#endif
+	min_log = log(bsdf_min*.5 + 1e-5);
 						/* output BSDF rep. */
 	for (n = 0; (n < 6) & (2*n+3 < argc); n++) {
-		dir[2] = sin((M_PI/180.)*atof(argv[2*n+2]));
-		dir[0] = dir[2] * cos((M_PI/180.)*atof(argv[2*n+3]));
-		dir[1] = dir[2] * sin((M_PI/180.)*atof(argv[2*n+3]));
-		dir[2] = input_orient * sqrt(1. - dir[2]*dir[2]);
-		fprintf(stderr, "Computing DSF for incident direction (%.1f,%.1f)\n",
-				get_theta180(dir), get_phi360(dir));
-		rbf = advect_rbf(dir, 15000);
-		if (rbf == NULL)
-			fputs("NULL RBF\n", stderr);
+		double	theta = atof(argv[2*n+2]);
+		if (inpXML) {
+			input_orient = (theta <= 90.) ? 1 : -1;
+			output_orient = doTrans ? -input_orient : input_orient;
+		}
+		idir[2] = sin((M_PI/180.)*theta);
+		idir[0] = idir[2] * cos((M_PI/180.)*atof(argv[2*n+3]));
+		idir[1] = idir[2] * sin((M_PI/180.)*atof(argv[2*n+3]));
+		idir[2] = input_orient * sqrt(1. - idir[2]*idir[2]);
+#ifdef DEBUG
+		fprintf(stderr, "Computing BSDF for incident direction (%.1f,%.1f)\n",
+				get_theta180(idir), get_phi360(idir));
+#endif
+		if (!inpXML)
+			rbf = advect_rbf(idir, 15000);
+#ifdef DEBUG
+		if (inpXML)
+			fprintf(stderr, "Hemispherical %s: %.3f\n",
+				(output_orient > 0 ? "reflection" : "transmission"),
+				SDdirectHemi(idir, SDsampSp|SDsampDf |
+						(output_orient > 0 ?
+						 SDsampR : SDsampT), &myBSDF));
+		else if (rbf == NULL)
+			fputs("Empty RBF\n", stderr);
 		else
-			fprintf(stderr, "Hemispherical reflectance: %.3f\n", rbf->vtotal);
+			fprintf(stderr, "Hemispherical %s: %.3f\n",
+				(output_orient > 0 ? "reflection" : "transmission"),
+				rbf->vtotal);
+#endif
 		printf("void trans tmat\n0\n0\n7 %f %f %f .04 .04 .9 1\n",
 				colarr[n][0], colarr[n][1], colarr[n][2]);
 		if (showPeaks && rbf != NULL) {
 			printf("void plastic pmat\n0\n0\n5 %f %f %f .04 .08\n",
 				1.-colarr[n][0], 1.-colarr[n][1], 1.-colarr[n][2]);
 			for (i = 0; i < rbf->nrbf; i++) {
-				ovec_from_pos(dir, rbf->rbfa[i].gx, rbf->rbfa[i].gy);
-				bsdf = eval_rbfrep(rbf, dir) / (output_orient*dir[2]);
-				bsdf = log(bsdf) - min_log;
+				ovec_from_pos(odir, rbf->rbfa[i].gx, rbf->rbfa[i].gy);
+				bsdf = eval_rbfrep(rbf, odir) / (output_orient*odir[2]);
+				bsdf = log(bsdf + 1e-5) - min_log;
 				printf("pmat sphere p%d\n0\n0\n4 %f %f %f %f\n",
-					i+1, dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf,
+					i+1, odir[0]*bsdf, odir[1]*bsdf, odir[2]*bsdf,
 						.007*bsdf);
 			}
 		}
@@ -89,11 +146,19 @@ main(int argc, char *argv[])
 		}
 		for (i = 0; i < GRIDRES; i++)
 		    for (j = 0; j < GRIDRES; j++) {
-			ovec_from_pos(dir, i, j);
-			bsdf = eval_rbfrep(rbf, dir) / (output_orient*dir[2]);
-			bsdf = log(bsdf) - min_log;
+			ovec_from_pos(odir, i, j);
+			if (inpXML) {
+				SDValue	sval;
+				if (SDreportError(SDevalBSDF(&sval, odir,
+							idir, &myBSDF), stderr))
+					return(1);
+				bsdf = sval.cieY;
+			} else
+				bsdf = eval_rbfrep(rbf, odir) /
+						(output_orient*odir[2]);
+			bsdf = log(bsdf + 1e-5) - min_log;
 			fprintf(fp, "%.8e %.8e %.8e\n",
-					dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
+					odir[0]*bsdf, odir[1]*bsdf, odir[2]*bsdf);
 		    }
 		if (rbf != NULL)
 			free(rbf);