--- ray/src/cv/bsdf2klems.c	2013/04/23 23:19:09	2.3
+++ ray/src/cv/bsdf2klems.c	2013/08/02 20:56:19	2.8
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf2klems.c,v 2.3 2013/04/23 23:19:09 greg Exp $";
+static const char RCSid[] = "$Id: bsdf2klems.c,v 2.8 2013/08/02 20:56:19 greg Exp $";
 #endif
 /*
  * Load measured BSDF interpolant and write out as XML file with Klems matrix.
@@ -17,6 +17,8 @@ static const char RCSid[] = "$Id: bsdf2klems.c,v 2.3 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 */
@@ -152,6 +154,7 @@ static void
 eval_bsdf(const char *fname)
 {
 	ANGLE_BASIS	*abp = get_basis(kbasis);
+	float		*trans_mtx = NULL;
 	SDData		bsd;
 	SDError		ec;
 	FVECT		vin, vout;
@@ -172,7 +175,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;
@@ -193,7 +196,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;
@@ -207,6 +210,9 @@ eval_bsdf(const char *fname)
 	}
 						/* front transmission */
 	if (bsd.tf != NULL || bsd.tLamb.cieY > .002) {
+	    if (bsd.tb == NULL)
+		trans_mtx = (float *)malloc(sizeof(float) *
+						abp->nangles*abp->nangles);
 	    input_orient = 1; output_orient = -1;
 	    data_prologue();
 	    for (j = 0; j < abp->nangles; j++) {
@@ -214,40 +220,63 @@ 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;
 			sum += sv.cieY;
 		    }
 		    printf("\t%.3e\n", sum/npsamps);
+		    if (trans_mtx != NULL)
+			trans_mtx[j*abp->nangles + i] = sum/npsamps;
 		}
 		putchar('\n');			/* extra space between rows */
 	    }
 	    data_epilogue();
 	}
 						/* back transmission */
-	if (bsd.tb != NULL) {
+	if (bsd.tb != NULL || trans_mtx != NULL) {
+	    if (bsd.tf == NULL)
+		trans_mtx = (float *)malloc(sizeof(float) *
+						abp->nangles*abp->nangles);
 	    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 */
-		    for (n = npsamps; n-- > 0; ) {
-			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
-			bi_getvec(vin, i+(n+frandom())/npsamps, abp);
-			ec = SDevalBSDF(&sv, vout, vin, &bsd);
-			if (ec != SDEnone)
+	        for (i = 0; i < abp->nangles; i++)
+		    if (bsd.tb != NULL) {	/* use tb if we have it */
+			sum = 0;		/* average over patches */
+			for (n = npsamps; n-- > 0; ) {
+			    fo_getvec(vout, j+(n+frandom())/npsamps, abp);
+			    bi_getvec(vin, i+urand(n), abp);
+			    ec = SDevalBSDF(&sv, vout, vin, &bsd);
+			    if (ec != SDEnone)
 				goto err;
-			sum += sv.cieY;
+			    sum += sv.cieY;
+			}
+			printf("\t%.3e\n", sum/npsamps);
+			if (trans_mtx != NULL)
+			    trans_mtx[i*abp->nangles + j] = sum/npsamps;
+		    } else {			/* else transpose tf */
+			printf("\t%.3e\n", trans_mtx[i*abp->nangles + j]);
 		    }
-		    printf("\t%.3e\n", sum/npsamps);
-		}
 		putchar('\n');			/* extra space between rows */
 	    }
 	    data_epilogue();
 	}
+						/* derived front transmission */
+	if (bsd.tf == NULL && trans_mtx != NULL) {
+	    input_orient = 1; output_orient = -1;
+	    data_prologue();
+	    for (j = 0; j < abp->nangles; j++) {
+	        for (i = 0; i < abp->nangles; i++)
+		    printf("\t%.3e\n", trans_mtx[j*abp->nangles + i]);
+		putchar('\n');			/* extra space between rows */
+	    }
+	    data_epilogue();
+	}
 	SDfreeBSDF(&bsd);			/* all done */
+	if (trans_mtx != NULL)
+		free(trans_mtx);
 	return;
 err:
 	SDreportError(ec, stderr);
@@ -259,10 +288,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++) {
@@ -274,10 +305,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);
@@ -291,7 +328,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;
@@ -324,6 +360,8 @@ eval_rbf(void)
 		}
 		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++) {
@@ -332,7 +370,6 @@ eval_rbf(void)
 	    putchar('\n');
 	}
 	data_epilogue();			/* finish output */
-#undef MAXPATCHES
 }
 
 /* Read in BSDF and interpolate as Klems matrix representation */
@@ -383,8 +420,11 @@ 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;
 		}
+		++eclock;
 		xml_header(argc, argv);			/* start XML output */
 		xml_prologue(NULL);
 		if (dofwd) {