src/util/rttree_reduce.c

#ifndef lint
static const char RCSid[] = "$Id: rttree_reduce.c,v 2.9 2014/03/17 01:59:48 greg Exp $";
#endif
/*
 *  A utility called by genBSDF.pl to reduce tensor tree samples and output
 *  in a standard format as required by XML specification for variable
 *  resolution BSDF data.  We are not meant to be run by the user directly.
 */

#include "rtio.h"
#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 = 95.;          /* target culling percentile */

#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)]

/* Tensor tree node */
typedef struct ttree_s {
        float           vmin, vmax;     /* value extrema */
        float           vavg;           /* average */
        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)
{
        pn->kid = (TNODE *)calloc(1<<ttrank, sizeof(TNODE));
        if (pn->kid == NULL)
                error(SYSTEM, "out of memory in new_kids");
}

/* Free children for this node */
static void
free_kids(TNODE *pn)
{
        int     i;

        if (pn->kid == NULL)
                return;
        for (i = 1<<ttrank; i--; )
                free_kids(pn->kid+i);
        free(pn->kid);
        pn->kid = NULL;
}

/* Build a tensor tree starting from the given hypercube */
static void
build_tree(TNODE *tp, const int bmin[], int l2s)
{
        int     bkmin[4];
        int     i, j;

        tp->vmin = 1e20;
        tp->vmax = 0;
        tp->vavg = 0;
        if (l2s <= 1) {                 /* reached upper leaves */
                for (i = 1<<ttrank; i--; ) {
                        float   val;
                        for (j = ttrank; j--; )
                                bkmin[j] = bmin[j] + (i>>j & 1);
                        val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
                                : dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
                        if (val < tp->vmin)
                                tp->vmin = val;
                        if (val > tp->vmax)
                                tp->vmax = val;
                        tp->vavg += val;
                }
                tp->vavg /= (float)(1<<ttrank);
                                        /* record stats */
                i = (HISTLEN/HISTMAX) * var_measure(tp);
                if (i >= HISTLEN) i = HISTLEN-1;
                ++histo[i];
                return;
        }
        --l2s;                          /* else still branching */
        new_kids(tp);                   /* grow recursively */
        for (i = 1<<ttrank; i--; ) {
                for (j = ttrank; j--; )
                        bkmin[j] = bmin[j] + ((i>>j & 1)<<l2s);
                build_tree(tp->kid+i, bkmin, l2s);
                if (tp->kid[i].vmin < tp->vmin)
                        tp->vmin = tp->kid[i].vmin;
                if (tp->kid[i].vmax > tp->vmax)
                        tp->vmax = tp->kid[i].vmax;
                tp->vavg += tp->kid[i].vavg;
        }
        tp->vavg /= (float)(1<<ttrank);
}

/* Set our trimming threshold */
static void
set_threshold(void)
{
        int     hsum = 0;
        int     i;

        for (i = HISTLEN; i--; )
                hsum += histo[i];
        hsum = pctcull*.01 * (double)hsum;
        for (i = 0; hsum > 0; i++)
                hsum -= histo[i];
        tthresh = (HISTMAX/HISTLEN) * i;
}

/* Trim our tree according to the current threshold */
static void
trim_tree(TNODE *tp)
{
        if (tp->kid == NULL)
                return;
        if (above_threshold(tp)) {      /* keeping branches? */
                int     i = 1<<ttrank;
                while (i--)
                        trim_tree(tp->kid+i);
                return;
        }
        free_kids(tp);                  /* else trim at this point */
}

/* Print a tensor tree from the given hypercube */
static void
print_tree(const TNODE *tp, const int bmin[], int l2s)
{
        int     bkmin[4];
        int     i, j;

        for (j = log2g-l2s; j--; )      /* indent based on branch level */
                fputs("    ", stdout);
        fputc('{', stdout);             /* special case for upper leaves */
        if (l2s <= 1 && above_threshold(tp)) {
                for (i = 0; i < 1<<ttrank; i++) {
                        float   val;
                        for (j = ttrank; j--; )
                                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);
                }
                fputs(" }\n", stdout);
                return;
        }
        if (tp->kid == NULL) {          /* trimmed limb */
                printf(" %.6e }\n", tp->vavg);
                return;
        }
        --l2s;                          /* else still branching */
        fputc('\n', stdout);
        for (i = 0; i < 1<<ttrank; i++) {
                for (j = ttrank; j--; )
                        bkmin[j] = bmin[j] + ((i>>j & 1)<<l2s);
                print_tree(tp->kid+i, bkmin, l2s);
        }
        ++l2s;
        for (j = log2g-l2s; j--; )
                fputs("    ", stdout);
        fputs("}\n", stdout);
}

/* Read a row of data in ASCII */
static int
read_ascii(float *rowp, int n)
{
        int     n2go;

        if ((rowp == NULL) | (n <= 0))
                return(0);
        for (n2go = n; n2go; n2go--)
                if (scanf("%f", rowp++) != 1)
                        break;
        if (n2go)
                error(USER, "unexpected EOD on ascii input");
        return(n-n2go);
}

/* Read a row of float data */
static int
read_float(float *rowp, int n)
{
        int     nread;

        if ((rowp == NULL) | (n <= 0))
                return(0);
        nread = fread(rowp, sizeof(float), n, stdin);
        if (nread != n)
                error(USER, "unexpected EOF on float input");
        return(nread);
}

/* Read a row of double data */
static int
read_double(float *rowp, int n)
{
        static double   *rowbuf = NULL;
        static int      rblen = 0;
        int             nread, i;

        if ((rowp == NULL) | (n <= 0)) {
                if (rblen) {
                        free(rowbuf);
                        rowbuf = NULL; rblen = 0;
                }
                return(0);
        }
        if (rblen < n) {
                rowbuf = (double *)realloc(rowbuf, sizeof(double)*(rblen=n));
                if (rowbuf == NULL)
                        error(SYSTEM, "out of memory in read_double");
        }
        nread = fread(rowbuf, sizeof(double), n, stdin);
        if (nread != n)
                error(USER, "unexpected EOF on double input");
        for (i = 0; i < nread; i++)
                *rowp++ = rowbuf[i];
        return(nread);
}

/* Truncate any negative values to zero */
static void
noneg(float *varr, int n)
{
        int     nnan = 0;

        while (n-- > 0) {
#ifdef isnan
                if (isnan(*varr)) {
                        *varr = 0;
                        ++nnan;
                } else
#endif
                if (*varr < 0) *varr = 0;
                ++varr;
        }
        if (nnan)
                fprintf(stderr, "Warning: BSDF data contains %d NaN values\n",
                                nnan);
}

/* Load data array, filling zeroes for rank 3 demi-tensor */
static void
load_data(void)
{
        int     (*readf)(float *, int) = NULL;
        
        switch (infmt) {
        case 'a':
                readf = &read_ascii;
                break;
        case 'f':
                readf = &read_float;
                break;
        case 'd':
                readf = &read_double;
                break;
        default:
                error(COMMAND, "unsupported input format");
                break;
        }
        datarr = (float *)calloc(1<<(log2g*ttrank), sizeof(float));
        if (datarr == NULL)
                error(SYSTEM, "out of memory in load_data");
        if (ttrank == 3) {
                int     ix, ox;
                for (ix = 0; ix < 1<<(log2g-1); ix++)
                        for (ox = 0; ox < 1<<log2g; ox++)
                                (*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),
                                                1<<log2g);
        }
        (*readf)(NULL, 0);      /* releases any buffers */
        if (infmt == 'a') {
                int     c;
                while ((c = getc(stdin)) != EOF) {
                        switch (c) {
                        case ' ':
                        case '\t':
                        case '\r':
                        case '\n':
                                continue;
                        }
                        error(WARNING, "data past end of expected input");
                        break;
                }
        } else if (getc(stdin) != EOF)
                error(WARNING, "binary data past end of expected input");

        noneg(datarr, 1<<(log2g*ttrank));       /* take precautions */
}

/* Enforce reciprocity by averaging data values */
static void
do_reciprocity(void)
{
        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;
                case 'r':
                        ttrank = atoi(argv[++i]);
                        if (ttrank != 3 && ttrank != 4)
                                goto userr;
                        break;
                case 'g':
                        log2g = atoi(argv[++i]);
                        if (log2g <= 1)
                                goto userr;
                        break;
                case 't':
                        pctcull = atof(argv[++i]);
                        if ((pctcull < 0) | (pctcull >= 100.))
                                goto userr;
                        break;
                case 'f':
                        infmt = argv[i][2];
                        if (!infmt || strchr("afd", infmt) == NULL)
                                goto userr;
                        break;
                default:
                        goto userr;
                }
        if (i < argc-1)
                goto userr;
                                        /* load input data */
        if (i == argc-1 && freopen(argv[i], "r", stdin) == NULL) {
                sprintf(errmsg, "cannot open input file '%s'", argv[i]);
                error(SYSTEM, errmsg);
        }
        if (infmt != 'a')
                SET_FILE_BINARY(stdin);
#ifdef getc_unlocked                    /* avoid lock/unlock overhead */
        flockfile(stdin);
#endif
        load_data();
        if (recipavg)
                do_reciprocity();
        if (doheader) {
                for (i = 0; i < argc; i++) {
                        fputs(argv[i], stdout);
                        fputc(i < argc-1 ? ' ' : '\n', stdout);
                }
                fputc('\n', stdout);
        }
        gtree.kid = NULL;               /* create our tree */
        bmin[0] = bmin[1] = bmin[2] = bmin[3] = 0;
        build_tree(&gtree, bmin, log2g);
                                        /* compute threshold & trim tree */
        set_threshold();
        trim_tree(&gtree);
                                        /* format to stdout */
        print_tree(&gtree, bmin, log2g);
        /* Clean up isn't necessary for main()...
        free_kids(&gtree);
        free(datarr);
        */
        return(0);
userr:
        fprintf(stderr, "Usage: %s [-h][-a][-f{a|f|d}][-r {3|4}][-g log2grid][-t trim%%] [input]\n",
                        argv[0]);
        return(1);
}
Revision:	2.10
Committed:	Tue Mar 18 01:33:07 2014 UTC (10 years, 1 month ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R2P2, rad4R2, rad4R2P1
Changes since 2.9:	+27 -4 lines
Log Message:	Added protection against negative and NaN values on input
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.10	static const char RCSid[] = "$Id: rttree_reduce.c,v 2.9 2014/03/17 01:59:48 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* A utility called by genBSDF.pl to reduce tensor tree samples and output
6			* in a standard format as required by XML specification for variable
7			* resolution BSDF data. We are not meant to be run by the user directly.
8			*/
9
10			#include "rtio.h"
11			#include "rterror.h"
12			#include "platform.h"
13			#include <stdlib.h>
14	greg	2.3	#include <math.h>
15	greg	2.1
16			float datarr; / our loaded BSDF data array */
17			int ttrank = 4; /* tensor tree rank */
18			int log2g = 4; /* log2 of grid resolution */
19			int infmt = 'a'; /* input format ('a','f','d') */
20	greg	2.4	double pctcull = 95.; /* target culling percentile */
21	greg	2.2
22	greg	2.1	#define dval3(ix,ox,oy) datarr[((((ix)<<log2g)+(ox))<<log2g)+(oy)]
23			#define dval4(ix,iy,ox,oy) datarr[((((((ix)<<log2g)+(iy))<<log2g)+(ox))<<log2g)+(oy)]
24
25			/* Tensor tree node */
26			typedef struct ttree_s {
27			float vmin, vmax; /* value extrema */
28			float vavg; /* average */
29			struct ttree_s kid; / 2^ttrank children */
30			} TNODE;
31
32	greg	2.3	#define HISTLEN 300 /* histogram resolution */
33			#define HISTMAX 10. /* maximum recorded measure in histogram */
34
35			int histo[HISTLEN]; /* histogram freq. of variance measure */
36
37			double tthresh; /* acceptance threshold (TBD) */
38
39			#define var_measure(tp) ( ((tp)->vmax - (tp)->vmin) / \
40			(sqrt((tp)->vavg) + .03) )
41			#define above_threshold(tp) (var_measure(tp) > tthresh)
42
43	greg	2.1	/* Allocate a new set of children for the given node (no checks) */
44			static void
45			new_kids(TNODE *pn)
46			{
47			pn->kid = (TNODE *)calloc(1<<ttrank, sizeof(TNODE));
48			if (pn->kid == NULL)
49			error(SYSTEM, "out of memory in new_kids");
50			}
51
52			/* Free children for this node */
53			static void
54			free_kids(TNODE *pn)
55			{
56			int i;
57
58			if (pn->kid == NULL)
59			return;
60			for (i = 1<<ttrank; i--; )
61			free_kids(pn->kid+i);
62			free(pn->kid);
63			pn->kid = NULL;
64			}
65
66			/* Build a tensor tree starting from the given hypercube */
67			static void
68			build_tree(TNODE *tp, const int bmin[], int l2s)
69			{
70			int bkmin[4];
71			int i, j;
72
73			tp->vmin = 1e20;
74			tp->vmax = 0;
75			tp->vavg = 0;
76			if (l2s <= 1) { /* reached upper leaves */
77			for (i = 1<<ttrank; i--; ) {
78			float val;
79			for (j = ttrank; j--; )
80			bkmin[j] = bmin[j] + (i>>j & 1);
81			val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
82			: dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
83			if (val < tp->vmin)
84			tp->vmin = val;
85			if (val > tp->vmax)
86			tp->vmax = val;
87			tp->vavg += val;
88			}
89			tp->vavg /= (float)(1<<ttrank);
90	greg	2.2	/* record stats */
91	greg	2.3	i = (HISTLEN/HISTMAX) * var_measure(tp);
92	greg	2.2	if (i >= HISTLEN) i = HISTLEN-1;
93			++histo[i];
94	greg	2.1	return;
95			}
96			--l2s; /* else still branching */
97			new_kids(tp); /* grow recursively */
98			for (i = 1<<ttrank; i--; ) {
99			for (j = ttrank; j--; )
100			bkmin[j] = bmin[j] + ((i>>j & 1)<<l2s);
101			build_tree(tp->kid+i, bkmin, l2s);
102			if (tp->kid[i].vmin < tp->vmin)
103			tp->vmin = tp->kid[i].vmin;
104			if (tp->kid[i].vmax > tp->vmax)
105			tp->vmax = tp->kid[i].vmax;
106			tp->vavg += tp->kid[i].vavg;
107			}
108			tp->vavg /= (float)(1<<ttrank);
109	greg	2.2	}
110
111			/* Set our trimming threshold */
112			static void
113	greg	2.10	set_threshold(void)
114	greg	2.2	{
115			int hsum = 0;
116			int i;
117
118			for (i = HISTLEN; i--; )
119			hsum += histo[i];
120			hsum = pctcull.01 (double)hsum;
121			for (i = 0; hsum > 0; i++)
122			hsum -= histo[i];
123			tthresh = (HISTMAX/HISTLEN) * i;
124			}
125
126			/* Trim our tree according to the current threshold */
127			static void
128			trim_tree(TNODE *tp)
129			{
130			if (tp->kid == NULL)
131			return;
132			if (above_threshold(tp)) { /* keeping branches? */
133			int i = 1<<ttrank;
134			while (i--)
135			trim_tree(tp->kid+i);
136			return;
137			}
138			free_kids(tp); /* else trim at this point */
139	greg	2.1	}
140
141			/* Print a tensor tree from the given hypercube */
142			static void
143			print_tree(const TNODE *tp, const int bmin[], int l2s)
144			{
145			int bkmin[4];
146			int i, j;
147
148			for (j = log2g-l2s; j--; ) /* indent based on branch level */
149			fputs(" ", stdout);
150			fputc('{', stdout); /* special case for upper leaves */
151			if (l2s <= 1 && above_threshold(tp)) {
152			for (i = 0; i < 1<<ttrank; i++) {
153			float val;
154			for (j = ttrank; j--; )
155	greg	2.5	bkmin[j] = bmin[j] + (i>>(ttrank-1-j) & 1);
156	greg	2.1	val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
157			: dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
158			printf(" %.4e", val);
159			}
160			fputs(" }\n", stdout);
161			return;
162			}
163			if (tp->kid == NULL) { /* trimmed limb */
164			printf(" %.6e }\n", tp->vavg);
165			return;
166			}
167			--l2s; /* else still branching */
168			fputc('\n', stdout);
169			for (i = 0; i < 1<<ttrank; i++) {
170			for (j = ttrank; j--; )
171			bkmin[j] = bmin[j] + ((i>>j & 1)<<l2s);
172			print_tree(tp->kid+i, bkmin, l2s);
173			}
174			++l2s;
175			for (j = log2g-l2s; j--; )
176			fputs(" ", stdout);
177			fputs("}\n", stdout);
178			}
179
180			/* Read a row of data in ASCII */
181			static int
182			read_ascii(float *rowp, int n)
183			{
184			int n2go;
185
186			if ((rowp == NULL) \| (n <= 0))
187			return(0);
188			for (n2go = n; n2go; n2go--)
189			if (scanf("%f", rowp++) != 1)
190			break;
191			if (n2go)
192			error(USER, "unexpected EOD on ascii input");
193			return(n-n2go);
194			}
195
196			/* Read a row of float data */
197			static int
198			read_float(float *rowp, int n)
199			{
200			int nread;
201
202			if ((rowp == NULL) \| (n <= 0))
203			return(0);
204			nread = fread(rowp, sizeof(float), n, stdin);
205			if (nread != n)
206			error(USER, "unexpected EOF on float input");
207			return(nread);
208			}
209
210			/* Read a row of double data */
211			static int
212			read_double(float *rowp, int n)
213			{
214			static double *rowbuf = NULL;
215			static int rblen = 0;
216			int nread, i;
217
218			if ((rowp == NULL) \| (n <= 0)) {
219			if (rblen) {
220			free(rowbuf);
221			rowbuf = NULL; rblen = 0;
222			}
223			return(0);
224			}
225			if (rblen < n) {
226			rowbuf = (double )realloc(rowbuf, sizeof(double)(rblen=n));
227			if (rowbuf == NULL)
228			error(SYSTEM, "out of memory in read_double");
229			}
230			nread = fread(rowbuf, sizeof(double), n, stdin);
231			if (nread != n)
232			error(USER, "unexpected EOF on double input");
233			for (i = 0; i < nread; i++)
234			*rowp++ = rowbuf[i];
235			return(nread);
236			}
237
238	greg	2.10	/* Truncate any negative values to zero */
239			static void
240			noneg(float *varr, int n)
241			{
242			int nnan = 0;
243
244			while (n-- > 0) {
245			#ifdef isnan
246			if (isnan(*varr)) {
247			*varr = 0;
248			++nnan;
249			} else
250			#endif
251			if (varr < 0) varr = 0;
252			++varr;
253			}
254			if (nnan)
255			fprintf(stderr, "Warning: BSDF data contains %d NaN values\n",
256			nnan);
257			}
258
259	greg	2.1	/* Load data array, filling zeroes for rank 3 demi-tensor */
260			static void
261	greg	2.10	load_data(void)
262	greg	2.1	{
263			int (readf)(float , int) = NULL;
264
265			switch (infmt) {
266			case 'a':
267			readf = &read_ascii;
268			break;
269			case 'f':
270			readf = &read_float;
271			break;
272			case 'd':
273			readf = &read_double;
274			break;
275			default:
276			error(COMMAND, "unsupported input format");
277			break;
278			}
279			datarr = (float )calloc(1<<(log2gttrank), sizeof(float));
280			if (datarr == NULL)
281			error(SYSTEM, "out of memory in load_data");
282			if (ttrank == 3) {
283			int ix, ox;
284	greg	2.4	for (ix = 0; ix < 1<<(log2g-1); ix++)
285	greg	2.1	for (ox = 0; ox < 1<<log2g; ox++)
286	greg	2.5	(*readf)(datarr+(((ix<<log2g)+ox)<<log2g),
287	greg	2.4	1<<log2g);
288	greg	2.1	} else /* ttrank == 4 */ {
289			int ix, iy, ox;
290			for (ix = 0; ix < 1<<log2g; ix++)
291			for (iy = 0; iy < 1<<log2g; iy++)
292			for (ox = 0; ox < 1<<log2g; ox++)
293			(*readf)(datarr +
294	greg	2.5	(((((ix<<log2g)+iy)<<log2g)+ox)<<log2g),
295	greg	2.1	1<<log2g);
296			}
297			(readf)(NULL, 0); / releases any buffers */
298			if (infmt == 'a') {
299			int c;
300			while ((c = getc(stdin)) != EOF) {
301			switch (c) {
302			case ' ':
303			case '\t':
304			case '\r':
305			case '\n':
306			continue;
307			}
308			error(WARNING, "data past end of expected input");
309			break;
310			}
311			} else if (getc(stdin) != EOF)
312			error(WARNING, "binary data past end of expected input");
313	greg	2.10
314			noneg(datarr, 1<<(log2gttrank)); / take precautions */
315	greg	2.1	}
316
317	greg	2.6	/* Enforce reciprocity by averaging data values */
318			static void
319	greg	2.10	do_reciprocity(void)
320	greg	2.6	{
321	greg	2.8	const int siz = 1<<log2g;
322			float v1p, v2p;
323	greg	2.6
324			if (ttrank == 3) {
325			int ix, ox, oy;
326	greg	2.8	for (ix = 0; ix < siz>>1; ix++)
327			for (ox = 0; ox < siz; ox++)
328			for (oy = 0; oy < siz>>1; oy++) {
329	greg	2.6	v1p = &dval3(ix,ox,oy);
330	greg	2.8	v2p = &dval3(ix,ox,siz-1-oy);
331	greg	2.6	v1p = v2p = .5f( v1p + *v2p );
332			}
333			} else /* ttrank == 4 */ {
334			int ix, iy, ox, oy;
335	greg	2.8	for (ix = 1; ix < siz; ix++)
336			for (iy = 1; iy < siz; iy++)
337	greg	2.6	for (ox = 0; ox < ix; ox++)
338			for (oy = 0; oy < iy; oy++) {
339	greg	2.8	v1p = &dval4(siz-1-ix,siz-1-iy,ox,oy);
340			v2p = &dval4(siz-1-ox,siz-1-oy,ix,iy);
341	greg	2.6	v1p = v2p = .5f( v1p + *v2p );
342			}
343			}
344			}
345
346	greg	2.1	/* Load BSDF array, coalesce uniform regions and format as tensor tree */
347			int
348			main(int argc, char *argv[])
349			{
350			int doheader = 1;
351	greg	2.6	int recipavg = 0;
352	greg	2.1	int bmin[4];
353			TNODE gtree;
354			int i;
355			/* get options and parameters */
356			for (i = 1; i < argc && argv[i][0] == '-'; i++)
357			switch (argv[i][1]) {
358	greg	2.6	case 'a':
359			recipavg = !recipavg;
360			break;
361	greg	2.1	case 'h':
362			doheader = !doheader;
363			break;
364			case 'r':
365			ttrank = atoi(argv[++i]);
366			if (ttrank != 3 && ttrank != 4)
367			goto userr;
368			break;
369			case 'g':
370			log2g = atoi(argv[++i]);
371			if (log2g <= 1)
372			goto userr;
373			break;
374			case 't':
375	greg	2.2	pctcull = atof(argv[++i]);
376			if ((pctcull < 0) \| (pctcull >= 100.))
377	greg	2.1	goto userr;
378			break;
379			case 'f':
380			infmt = argv[i][2];
381			if (!infmt \|\| strchr("afd", infmt) == NULL)
382			goto userr;
383			break;
384			default:
385			goto userr;
386			}
387			if (i < argc-1)
388			goto userr;
389			/* load input data */
390	greg	2.9	if (i == argc-1 && freopen(argv[i], "r", stdin) == NULL) {
391	greg	2.1	sprintf(errmsg, "cannot open input file '%s'", argv[i]);
392			error(SYSTEM, errmsg);
393			}
394			if (infmt != 'a')
395			SET_FILE_BINARY(stdin);
396	greg	2.9	#ifdef getc_unlocked /* avoid lock/unlock overhead */
397			flockfile(stdin);
398			#endif
399	greg	2.1	load_data();
400	greg	2.6	if (recipavg)
401			do_reciprocity();
402	greg	2.1	if (doheader) {
403			for (i = 0; i < argc; i++) {
404			fputs(argv[i], stdout);
405			fputc(i < argc-1 ? ' ' : '\n', stdout);
406			}
407			fputc('\n', stdout);
408			}
409			gtree.kid = NULL; /* create our tree */
410			bmin[0] = bmin[1] = bmin[2] = bmin[3] = 0;
411			build_tree(&gtree, bmin, log2g);
412	greg	2.2	/* compute threshold & trim tree */
413			set_threshold();
414			trim_tree(&gtree);
415	greg	2.1	/* format to stdout */
416			print_tree(&gtree, bmin, log2g);
417			/* Clean up isn't necessary for main()...
418			free_kids(&gtree);
419			free(datarr);
420			*/
421			return(0);
422			userr:
423	greg	2.7	fprintf(stderr, "Usage: %s [-h][-a][-f{a\|f\|d}][-r {3\|4}][-g log2grid][-t trim%%] [input]\n",
424	greg	2.1	argv[0]);
425			return(1);
426			}