src/util/rttree_reduce.c

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