--- ray/src/cv/bsdf2klems.c	2013/04/23 04:40:23	2.2
+++ ray/src/cv/bsdf2klems.c	2013/11/21 20:09:06	2.12
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf2klems.c,v 2.2 2013/04/23 04:40:23 greg Exp $";
+static const char RCSid[] = "$Id: bsdf2klems.c,v 2.12 2013/11/21 20:09:06 greg Exp $";
 #endif
 /*
  * Load measured BSDF interpolant and write out as XML file with Klems matrix.
@@ -17,12 +17,16 @@ static const char RCSid[] = "$Id: bsdf2klems.c,v 2.2 2
 #include "calcomp.h"
 #include "bsdfrep.h"
 #include "bsdf_m.h"
+				/* assumed maximum # Klems patches */
+#define MAXPATCHES	145
 				/* global argv[0] */
 char			*progname;
 				/* selected basis function name */
 static const char	*kbasis = "LBNL/Klems Full";
 				/* number of BSDF samples per patch */
 static int		npsamps = 256;
+				/* limit on number of RBF lobes */
+static int		lobe_lim = 15000;
 
 /* Return angle basis corresponding to the given name */
 ANGLE_BASIS *
@@ -36,17 +40,10 @@ get_basis(const char *bn)
 	return NULL;
 }
 
-/* Output XML prologue to stdout */
+/* Output XML header to stdout */
 static void
-xml_prologue(int ac, char *av[])
+xml_header(int ac, char *av[])
 {
-	ANGLE_BASIS	*abp = get_basis(kbasis);
-	int		i;
-
-	if (abp == NULL) {
-		fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis);
-		exit(1);
-	}
 	puts("<?xml version=\"1.0\" encoding=\"UTF-8\"?>");
 	puts("<WindowElement xmlns=\"http://windows.lbl.gov\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xsi:schemaLocation=\"http://windows.lbl.gov/BSDF-v1.4.xsd\">");
 	fputs("<!-- File produced by:", stdout);
@@ -55,15 +52,44 @@ xml_prologue(int ac, char *av[])
 		fputs(*av++, stdout);
 	}
 	puts(" -->");
+}
+
+/* Output XML prologue to stdout */
+static void
+xml_prologue(const SDData *sd)
+{
+	const char	*matn = (sd && sd->matn[0]) ? sd->matn :
+				bsdf_name[0] ? bsdf_name : "Unknown";
+	const char	*makr = (sd && sd->makr[0]) ? sd->makr :
+				bsdf_manuf[0] ? bsdf_manuf : "Unknown";
+	ANGLE_BASIS	*abp = get_basis(kbasis);
+	int		i;
+
+	if (abp == NULL) {
+		fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis);
+		exit(1);
+	}
 	puts("<WindowElementType>System</WindowElementType>");
 	puts("<FileType>BSDF</FileType>");
 	puts("<Optical>");
 	puts("<Layer>");
 	puts("\t<Material>");
-	puts("\t\t<Name>Name</Name>");
-	puts("\t\t<Manufacturer>Manufacturer</Manufacturer>");
+	printf("\t\t<Name>%s</Name>\n", matn);
+	printf("\t\t<Manufacturer>%s</Manufacturer>\n", makr);
+	if (sd && sd->dim[2] > .001) {
+		printf("\t\t<Thickness unit=\"meter\">%.3f</Thickness>\n", sd->dim[2]);
+		printf("\t\t<Width unit=\"meter\">%.3f</Width>\n", sd->dim[0]);
+		printf("\t\t<Height unit=\"meter\">%.3f</Height>\n", sd->dim[1]);
+	}
 	puts("\t\t<DeviceType>Other</DeviceType>");
 	puts("\t</Material>");
+	if (sd && sd->mgf != NULL) {
+		puts("\t<Geometry format=\"MGF\">");
+		puts("\t\t<MGFblock unit=\"meter\">");
+		fputs(sd->mgf, stdout);
+		puts("</MGFblock>");
+		puts("\t</Geometry>");
+	}
 	puts("\t<DataDefinition>");
 	puts("\t\t<IncidentDataStructure>Columns</IncidentDataStructure>");
 	puts("\t\t<AngleBasis>");
@@ -73,7 +99,7 @@ xml_prologue(int ac, char *av[])
 		printf("\t\t\t<Theta>%g</Theta>\n", i ?
 				.5*(abp->lat[i].tmin + abp->lat[i+1].tmin) :
 				.0 );
-		printf("\t\t\t<nPhis>%d</nPhis>", abp->lat[i].nphis);
+		printf("\t\t\t<nPhis>%d</nPhis>\n", abp->lat[i].nphis);
 		puts("\t\t\t<ThetaBounds>");
 		printf("\t\t\t\t<LowerTheta>%g</LowerTheta>\n", abp->lat[i].tmin);
 		printf("\t\t\t\t<UpperTheta>%g</UpperTheta>\n", abp->lat[i+1].tmin);
@@ -142,6 +168,7 @@ eval_bsdf(const char *fname)
 	SDclearBSDF(&bsd, fname);		/* load BSDF file */
 	if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
 		goto err;
+	xml_prologue(&bsd);			/* pass geometry */
 						/* front reflection */
 	if (bsd.rf != NULL || bsd.rLambFront.cieY > .002) {
 	    input_orient = 1; output_orient = 1;
@@ -151,7 +178,7 @@ eval_bsdf(const char *fname)
 		    sum = 0;			/* average over patches */
 		    for (n = npsamps; n-- > 0; ) {
 			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
-			fi_getvec(vin, i+(n+frandom())/npsamps, abp);
+			fi_getvec(vin, i+urand(n), abp);
 			ec = SDevalBSDF(&sv, vout, vin, &bsd);
 			if (ec != SDEnone)
 				goto err;
@@ -172,7 +199,7 @@ eval_bsdf(const char *fname)
 		    sum = 0;			/* average over patches */
 		    for (n = npsamps; n-- > 0; ) {
 			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
-			bi_getvec(vin, i+(n+frandom())/npsamps, abp);
+			bi_getvec(vin, i+urand(n), abp);
 			ec = SDevalBSDF(&sv, vout, vin, &bsd);
 			if (ec != SDEnone)
 				goto err;
@@ -193,7 +220,7 @@ eval_bsdf(const char *fname)
 		    sum = 0;			/* average over patches */
 		    for (n = npsamps; n-- > 0; ) {
 			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
-			fi_getvec(vin, i+(n+frandom())/npsamps, abp);
+			fi_getvec(vin, i+urand(n), abp);
 			ec = SDevalBSDF(&sv, vout, vin, &bsd);
 			if (ec != SDEnone)
 				goto err;
@@ -206,15 +233,15 @@ eval_bsdf(const char *fname)
 	    data_epilogue();
 	}
 						/* back transmission */
-	if (bsd.tb != NULL) {
+	if ((bsd.tb != NULL) | (bsd.tf != NULL)) {
 	    input_orient = -1; output_orient = 1;
 	    data_prologue();
 	    for (j = 0; j < abp->nangles; j++) {
 	        for (i = 0; i < abp->nangles; i++) {
-		    sum = 0;			/* average over patches */
+		    sum = 0;		/* average over patches */
 		    for (n = npsamps; n-- > 0; ) {
 			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
-			bi_getvec(vin, i+(n+frandom())/npsamps, abp);
+			bi_getvec(vin, i+urand(n), abp);
 			ec = SDevalBSDF(&sv, vout, vin, &bsd);
 			if (ec != SDEnone)
 				goto err;
@@ -238,10 +265,12 @@ static void
 eval_function(char *funame)
 {
 	ANGLE_BASIS	*abp = get_basis(kbasis);
+	int		assignD = (fundefined(funame) < 6);
 	double		iovec[6];
 	double		sum;
 	int		i, j, n;
 
+	initurand(npsamps);
 	data_prologue();			/* begin output */
 	for (j = 0; j < abp->nangles; j++) {	/* run through directions */
 	    for (i = 0; i < abp->nangles; i++) {
@@ -253,10 +282,16 @@ eval_function(char *funame)
 			bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
 
 		    if (input_orient > 0)
-			fi_getvec(iovec, i+(n+frandom())/npsamps, abp);
+			fi_getvec(iovec, i+urand(n), abp);
 		    else
-			bi_getvec(iovec, i+(n+frandom())/npsamps, abp);
+			bi_getvec(iovec, i+urand(n), abp);
 
+		    if (assignD) {
+			varset("Dx", '=', -iovec[3]);
+			varset("Dy", '=', -iovec[4]);
+			varset("Dz", '=', -iovec[5]);
+			++eclock;
+		    }
 		    sum += funvalue(funame, 6, iovec);
 		}
 		printf("\t%.3e\n", sum/npsamps);
@@ -270,7 +305,6 @@ eval_function(char *funame)
 static void
 eval_rbf(void)
 {
-#define MAXPATCHES	145
 	ANGLE_BASIS	*abp = get_basis(kbasis);
 	float		bsdfarr[MAXPATCHES*MAXPATCHES];
 	FVECT		vin, vout;
@@ -289,7 +323,7 @@ eval_rbf(void)
 	    else
 		bi_getvec(vin, i+.5*(i>0), abp);
 
-	    rbf = advect_rbf(vin);		/* compute radial basis func */
+	    rbf = advect_rbf(vin, lobe_lim);	/* compute radial basis func */
 
 	    for (j = 0; j < abp->nangles; j++) {
 	        sum = 0;			/* sample over exiting patch */
@@ -299,10 +333,13 @@ eval_rbf(void)
 		    else
 			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
 
-		    sum += eval_rbfrep(rbf, vout) / vout[2];
+		    sum += eval_rbfrep(rbf, vout);
 		}
-		bsdfarr[j*abp->nangles + i] = sum*output_orient/npsamps;
+		fo_getvec(vout, j+.5, abp);	/* use centered secant */
+		bsdfarr[j*abp->nangles + i] = sum / (npsamps*vout[2]);
 	    }
+	    if (rbf != NULL)
+		free(rbf);
 	}
 	n = 0;					/* write out our matrix */
 	for (j = 0; j < abp->nangles; j++) {
@@ -311,7 +348,6 @@ eval_rbf(void)
 	    putchar('\n');
 	}
 	data_epilogue();			/* finish output */
-#undef MAXPATCHES
 }
 
 /* Read in BSDF and interpolate as Klems matrix representation */
@@ -354,6 +390,9 @@ main(int argc, char *argv[])
 		case 'q':
 			kbasis = "LBNL/Klems Quarter";
 			break;
+		case 'l':
+			lobe_lim = atoi(argv[++i]);
+			break;
 		default:
 			goto userr;
 		}
@@ -362,9 +401,13 @@ main(int argc, char *argv[])
 			fprintf(stderr,
 	"%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
 					progname);
+			fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n",
+					progname);
 			goto userr;
 		}
-		xml_prologue(argc, argv);	/* start XML output */
+		++eclock;
+		xml_header(argc, argv);			/* start XML output */
+		xml_prologue(NULL);
 		if (dofwd) {
 			input_orient = -1;
 			output_orient = -1;
@@ -385,7 +428,7 @@ main(int argc, char *argv[])
 						/* XML input? */
 	if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
 				!strcasecmp(cp, ".xml")) {
-		xml_prologue(argc, argv);	/* start XML output */
+		xml_header(argc, argv);		/* start XML output */
 		eval_bsdf(argv[i]);		/* load & resample BSDF */
 		xml_epilogue();			/* finish XML output & exit */
 		return(0);
@@ -402,8 +445,10 @@ main(int argc, char *argv[])
 			if (!load_bsdf_rep(fpin))
 				return(1);
 			fclose(fpin);
-			if (!nbsdf++)		/* start XML on first dist. */
-				xml_prologue(argc, argv);
+			if (!nbsdf++) {		/* start XML on first dist. */
+				xml_header(argc, argv);
+				xml_prologue(NULL);
+			}
 			eval_rbf();
 		}
 		xml_epilogue();			/* finish XML output & exit */
@@ -412,17 +457,16 @@ main(int argc, char *argv[])
 	SET_FILE_BINARY(stdin);			/* load from stdin */
 	if (!load_bsdf_rep(stdin))
 		return(1);
-	xml_prologue(argc, argv);		/* start XML output */
+	xml_header(argc, argv);			/* start XML output */
+	xml_prologue(NULL);
 	eval_rbf();				/* resample dist. */
 	xml_epilogue();				/* finish XML output & exit */
 	return(0);
 userr:
 	fprintf(stderr,
-	"Usage: %s [-n spp][-h|-q][bsdf.sir ..] > bsdf.xml\n",
-				progname);
+	"Usage: %s [-n spp][-h|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
 	fprintf(stderr,
-	"   or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n",
-				progname);
+	"   or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
 	fprintf(stderr,
 	"   or: %s [-n spp][-h|-q][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
 				progname);