--- ray/src/util/rttree_reduce.c	2011/05/31 20:50:26	2.2
+++ ray/src/util/rttree_reduce.c	2013/10/03 03:36:11	2.8
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: rttree_reduce.c,v 2.2 2011/05/31 20:50:26 greg Exp $";
+static const char RCSid[] = "$Id: rttree_reduce.c,v 2.8 2013/10/03 03:36:11 greg Exp $";
 #endif
 /*
  *  A utility called by genBSDF.pl to reduce tensor tree samples and output
@@ -11,25 +11,17 @@ static const char RCSid[] = "$Id: rttree_reduce.c,v 2.
 #include "rterror.h"
 #include "platform.h"
 #include <stdlib.h>
+#include <math.h>
 
 float	*datarr;		/* our loaded BSDF data array */
 int	ttrank = 4;		/* tensor tree rank */
 int	log2g = 4;		/* log2 of grid resolution */
 int	infmt = 'a';		/* input format ('a','f','d') */
-double	pctcull = 99.;		/* target culling percentile */
+double	pctcull = 95.;		/* target culling percentile */
 
-#define HISTLEN		300	/* histogram resolution */
-#define HISTMAX		10.	/* maximum recorded value in histogram */
-
-int	histo[HISTLEN];		/* histogram freq. of max-min BSDF */
-
-double	tthresh;		/* acceptance threshold (TBD) */
-
 #define dval3(ix,ox,oy)		datarr[((((ix)<<log2g)+(ox))<<log2g)+(oy)]
 #define dval4(ix,iy,ox,oy)	datarr[((((((ix)<<log2g)+(iy))<<log2g)+(ox))<<log2g)+(oy)]
 
-#define above_threshold(tp)	((tp)->vmax - (tp)->vmin > tthresh)
-
 /* Tensor tree node */
 typedef struct ttree_s {
 	float		vmin, vmax;	/* value extrema */
@@ -37,6 +29,17 @@ typedef struct ttree_s {
 	struct ttree_s	*kid;		/* 2^ttrank children */
 } TNODE;
 
+#define HISTLEN		300	/* histogram resolution */
+#define HISTMAX		10.	/* maximum recorded measure in histogram */
+
+int	histo[HISTLEN];		/* histogram freq. of variance measure */
+
+double	tthresh;		/* acceptance threshold (TBD) */
+
+#define var_measure(tp)		( ((tp)->vmax - (tp)->vmin) / \
+					(sqrt((tp)->vavg) + .03) )
+#define above_threshold(tp)	(var_measure(tp) > tthresh)
+
 /* Allocate a new set of children for the given node (no checks) */
 static void
 new_kids(TNODE *pn)
@@ -85,7 +88,7 @@ build_tree(TNODE *tp, const int bmin[], int l2s)
 		}
 		tp->vavg /= (float)(1<<ttrank);
 					/* record stats */
-		i = (HISTLEN/HISTMAX) * (tp->vmax - tp->vmin);
+		i = (HISTLEN/HISTMAX) * var_measure(tp);
 		if (i >= HISTLEN) i = HISTLEN-1;
 		++histo[i];
 		return;
@@ -149,7 +152,7 @@ print_tree(const TNODE *tp, const int bmin[], int l2s)
 		for (i = 0; i < 1<<ttrank; i++) {
 			float	val;
 			for (j = ttrank; j--; )
-				bkmin[j] = bmin[j] + (i>>j & 1);
+				bkmin[j] = bmin[j] + (i>>(ttrank-1-j) & 1);
 			val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
 				: dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
 			printf(" %.4e", val);
@@ -257,17 +260,17 @@ load_data()
 		error(SYSTEM, "out of memory in load_data");
 	if (ttrank == 3) {
 		int	ix, ox;
-		for (ix = 0; ix < 1<<log2g; ix++)
+		for (ix = 0; ix < 1<<(log2g-1); ix++)
 			for (ox = 0; ox < 1<<log2g; ox++)
-				(*readf)(datarr+((((ix)<<log2g)+(ox))<<log2g),
-						1<<(log2g-1));
+				(*readf)(datarr+(((ix<<log2g)+ox)<<log2g),
+						1<<log2g);
 	} else /* ttrank == 4 */ {
 		int	ix, iy, ox;
 		for (ix = 0; ix < 1<<log2g; ix++)
 		    for (iy = 0; iy < 1<<log2g; iy++)
 			for (ox = 0; ox < 1<<log2g; ox++)
 				(*readf)(datarr +
-				((((((ix)<<log2g)+(iy))<<log2g)+(ox))<<log2g),
+				(((((ix<<log2g)+iy)<<log2g)+ox)<<log2g),
 						1<<log2g);
 	}
 	(*readf)(NULL, 0);	/* releases any buffers */
@@ -288,17 +291,50 @@ load_data()
 		error(WARNING, "binary data past end of expected input");
 }
 
+/* Enforce reciprocity by averaging data values */
+static void
+do_reciprocity()
+{
+	const int	siz = 1<<log2g;
+	float		*v1p, *v2p;
+
+	if (ttrank == 3) {
+		int	ix, ox, oy;
+		for (ix = 0; ix < siz>>1; ix++)
+		    for (ox = 0; ox < siz; ox++)
+		        for (oy = 0; oy < siz>>1; oy++) {
+				v1p = &dval3(ix,ox,oy);
+				v2p = &dval3(ix,ox,siz-1-oy);
+				*v1p = *v2p = .5f*( *v1p + *v2p );
+			}
+	} else /* ttrank == 4 */ {
+		int	ix, iy, ox, oy;
+		for (ix = 1; ix < siz; ix++)
+		    for (iy = 1; iy < siz; iy++)
+			for (ox = 0; ox < ix; ox++)
+			    for (oy = 0; oy < iy; oy++) {
+				v1p = &dval4(siz-1-ix,siz-1-iy,ox,oy);
+				v2p = &dval4(siz-1-ox,siz-1-oy,ix,iy);
+				*v1p = *v2p = .5f*( *v1p + *v2p );
+			    }
+	}
+}
+
 /* Load BSDF array, coalesce uniform regions and format as tensor tree */
 int
 main(int argc, char *argv[])
 {
 	int	doheader = 1;
+	int	recipavg = 0;
 	int	bmin[4];
 	TNODE	gtree;
 	int	i;
 					/* get options and parameters */
 	for (i = 1; i < argc && argv[i][0] == '-'; i++)
 		switch (argv[i][1]) {
+		case 'a':
+			recipavg = !recipavg;
+			break;
 		case 'h':
 			doheader = !doheader;
 			break;
@@ -335,6 +371,8 @@ main(int argc, char *argv[])
 	if (infmt != 'a')
 		SET_FILE_BINARY(stdin);
 	load_data();
+	if (recipavg)
+		do_reciprocity();
 	if (doheader) {
 		for (i = 0; i < argc; i++) {
 			fputs(argv[i], stdout);
@@ -356,7 +394,7 @@ main(int argc, char *argv[])
 	*/
 	return(0);
 userr:
-	fprintf(stderr, "Usage: %s [-h][-f{a|f|d}][-r {3|4}][-g log2grid][-t trim%%] [input]\n",
+	fprintf(stderr, "Usage: %s [-h][-a][-f{a|f|d}][-r {3|4}][-g log2grid][-t trim%%] [input]\n",
 			argv[0]);
 	return(1);
 }