--- ray/src/cv/bsdf2klems.c	2013/04/21 23:01:14	2.1
+++ ray/src/cv/bsdf2klems.c	2013/11/09 05:47:49	2.11
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf2klems.c,v 2.1 2013/04/21 23:01:14 greg Exp $";
+static const char RCSid[] = "$Id: bsdf2klems.c,v 2.11 2013/11/09 05:47:49 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.1 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);
@@ -127,7 +153,7 @@ xml_epilogue(void)
 	puts("</WindowElement>");
 }
 
-/* Load and resample XML BSDF description */
+/* Load and resample XML BSDF description using Klems basis */
 static void
 eval_bsdf(const char *fname)
 {
@@ -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;
@@ -233,15 +260,17 @@ err:
 	exit(1);
 }
 
-/* Interpolate and output a BSDF function */
+/* Interpolate and output a BSDF function using Klems basis */
 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, j+(n+frandom())/npsamps, abp);
+			fi_getvec(iovec, i+urand(n), abp);
 		    else
-			bi_getvec(iovec, j+(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);
@@ -271,12 +306,47 @@ static void
 eval_rbf(void)
 {
 	ANGLE_BASIS	*abp = get_basis(kbasis);
-	double		iovec[6];
+	float		bsdfarr[MAXPATCHES*MAXPATCHES];
+	FVECT		vin, vout;
+	RBFNODE		*rbf;
 	double		sum;
 	int		i, j, n;
+						/* sanity check */
+	if (abp->nangles > MAXPATCHES) {
+		fprintf(stderr, "%s: too many patches!\n", progname);
+		exit(1);
+	}
+	data_prologue();			/* begin output */
+	for (i = 0; i < abp->nangles; i++) {
+	    if (input_orient > 0)		/* use incident patch center */
+		fi_getvec(vin, i+.5*(i>0), abp);
+	    else
+		bi_getvec(vin, i+.5*(i>0), abp);
 
-fprintf(stder, "%s: RBF evaluation currently unimplemented\n", progname);
-exit(1);
+	    rbf = advect_rbf(vin, lobe_lim);	/* compute radial basis func */
+
+	    for (j = 0; j < abp->nangles; j++) {
+	        sum = 0;			/* sample over exiting patch */
+		for (n = npsamps; n--; ) {
+		    if (output_orient > 0)
+			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
+		    else
+			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
+
+		    sum += eval_rbfrep(rbf, vout) / vout[2];
+		}
+		bsdfarr[j*abp->nangles + i] = sum*output_orient/npsamps;
+	    }
+	    if (rbf != NULL)
+		free(rbf);
+	}
+	n = 0;					/* write out our matrix */
+	for (j = 0; j < abp->nangles; j++) {
+	    for (i = 0; i < abp->nangles; i++)
+		printf("\t%.3e\n", bsdfarr[n++]);
+	    putchar('\n');
+	}
+	data_epilogue();			/* finish output */
 }
 
 /* Read in BSDF and interpolate as Klems matrix representation */
@@ -292,7 +362,7 @@ main(int argc, char *argv[])
 	esupport &= ~(E_INCHAN|E_OUTCHAN);
 	scompile("PI:3.14159265358979323846", NULL, 0);
 	biggerlib();
-	for (i = 1; i < argc-1 && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
+	for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
 		switch (argv[i][1]) {		/* get options */
 		case 'n':
 			npsamps = atoi(argv[++i]);
@@ -319,6 +389,9 @@ main(int argc, char *argv[])
 		case 'q':
 			kbasis = "LBNL/Klems Quarter";
 			break;
+		case 'l':
+			lobe_lim = atoi(argv[++i]);
+			break;
 		default:
 			goto userr;
 		}
@@ -327,9 +400,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;
@@ -347,9 +424,10 @@ main(int argc, char *argv[])
 		xml_epilogue();			/* finish XML output & exit */
 		return(0);
 	}
-	if (i == argc-1 && (cp = strstr(argv[i], ".xml")) != NULL &&
-			strlen(cp) == 4) {	/* XML input? */
-		xml_prologue(argc, argv);	/* start XML output */
+						/* XML input? */
+	if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
+				!strcasecmp(cp, ".xml")) {
+		xml_header(argc, argv);		/* start XML output */
 		eval_bsdf(argv[i]);		/* load & resample BSDF */
 		xml_epilogue();			/* finish XML output & exit */
 		return(0);
@@ -366,8 +444,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 */
@@ -376,17 +456,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);