src/util/rttree_reduce.c

#ifndef lint
static const char RCSid[] = "$Id: rttree_reduce.c,v 2.18 2019/03/02 17:05:45 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[((size_t)(((ix)<<log2g)+(ox))<<log2g)+(oy)]
#define dval4(ix,iy,ox,oy)      datarr[((((size_t)(((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         4.      /* maximum recorded measure in histogram */

int     histo[HISTLEN];         /* histogram freq. of variance measure */

double  tthresh;                /* acceptance threshold (TBD) */

#define var_measure(tp)         sqrt( ((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((size_t)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 = getbinary(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) {
                if (rblen) free(rowbuf);
                rowbuf = (double *)malloc(sizeof(double)*(rblen=n));
                if (rowbuf == NULL)
                        error(SYSTEM, "out of memory in read_double");
        }
        nread = getbinary(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((size_t)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)(&dval3(ix,ox,0), 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)(&dval4(ix,iy,ox,0), 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 = 0; ix < siz; ix++)
                    for (iy = 0; iy < siz; iy++) {
                        int     cnt = ix*siz + iy;
                        for (ox = 0; cnt > 0; ox++)
                            for (oy = 0; oy < siz; 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 );
                                if (--cnt <= 0)
                                        break;
                            }
                    }
        }
}

/* 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.19
Committed:	Mon May 11 20:26:12 2020 UTC (5 years, 5 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.18:	+3 -3 lines
Log Message:	Fixed integer overflow for very large trees
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.19	static const char RCSid[] = "$Id: rttree_reduce.c,v 2.18 2019/03/02 17:05:45 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.19	#define dval3(ix,ox,oy) datarr[((size_t)(((ix)<<log2g)+(ox))<<log2g)+(oy)]
23			#define dval4(ix,iy,ox,oy) datarr[((((size_t)(((ix)<<log2g)+(iy))<<log2g)+(ox))<<log2g)+(oy)]
24	greg	2.1
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	greg	2.18	#define HISTMAX 4. /* maximum recorded measure in histogram */
34	greg	2.3
35			int histo[HISTLEN]; /* histogram freq. of variance measure */
36
37			double tthresh; /* acceptance threshold (TBD) */
38
39	greg	2.18	#define var_measure(tp) sqrt( ((tp)->vmax - (tp)->vmin) / \
40	greg	2.3	(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	greg	2.14	pn->kid = (TNODE *)calloc((size_t)1<<ttrank, sizeof(TNODE));
48	greg	2.1	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	greg	2.12	printf((0.001<=val)&(val<10.) ? " %.7f" : " %.3e", val);
159	greg	2.1	}
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	greg	2.16	nread = getbinary(rowp, sizeof(float), n, stdin);
205	greg	2.1	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	greg	2.17	if (rblen) free(rowbuf);
227			rowbuf = (double )malloc(sizeof(double)(rblen=n));
228	greg	2.1	if (rowbuf == NULL)
229			error(SYSTEM, "out of memory in read_double");
230			}
231	greg	2.16	nread = getbinary(rowbuf, sizeof(double), n, stdin);
232	greg	2.1	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.15	datarr = (float )calloc((size_t)1<<(log2gttrank), sizeof(float));
281	greg	2.1	if (datarr == NULL)
282			error(SYSTEM, "out of memory in load_data");
283			if (ttrank == 3) {
284			int ix, ox;
285	greg	2.4	for (ix = 0; ix < 1<<(log2g-1); ix++)
286	greg	2.1	for (ox = 0; ox < 1<<log2g; ox++)
287	greg	2.17	(*readf)(&dval3(ix,ox,0), 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	greg	2.17	(*readf)(&dval4(ix,iy,ox,0), 1<<log2g);
294	greg	2.1	}
295			(readf)(NULL, 0); / releases any buffers */
296			if (infmt == 'a') {
297			int c;
298			while ((c = getc(stdin)) != EOF) {
299			switch (c) {
300			case ' ':
301			case '\t':
302			case '\r':
303			case '\n':
304			continue;
305			}
306			error(WARNING, "data past end of expected input");
307			break;
308			}
309			} else if (getc(stdin) != EOF)
310			error(WARNING, "binary data past end of expected input");
311	greg	2.10
312			noneg(datarr, 1<<(log2gttrank)); / take precautions */
313	greg	2.1	}
314
315	greg	2.6	/* Enforce reciprocity by averaging data values */
316			static void
317	greg	2.10	do_reciprocity(void)
318	greg	2.6	{
319	greg	2.8	const int siz = 1<<log2g;
320			float v1p, v2p;
321	greg	2.6
322			if (ttrank == 3) {
323			int ix, ox, oy;
324	greg	2.8	for (ix = 0; ix < siz>>1; ix++)
325			for (ox = 0; ox < siz; ox++)
326			for (oy = 0; oy < siz>>1; oy++) {
327	greg	2.6	v1p = &dval3(ix,ox,oy);
328	greg	2.8	v2p = &dval3(ix,ox,siz-1-oy);
329	greg	2.6	v1p = v2p = .5f( v1p + *v2p );
330			}
331			} else /* ttrank == 4 */ {
332			int ix, iy, ox, oy;
333	greg	2.17	for (ix = 0; ix < siz; ix++)
334			for (iy = 0; iy < siz; iy++) {
335			int cnt = ix*siz + iy;
336			for (ox = 0; cnt > 0; ox++)
337			for (oy = 0; oy < siz; oy++) {
338	greg	2.8	v1p = &dval4(siz-1-ix,siz-1-iy,ox,oy);
339			v2p = &dval4(siz-1-ox,siz-1-oy,ix,iy);
340	greg	2.6	v1p = v2p = .5f( v1p + *v2p );
341	greg	2.17	if (--cnt <= 0)
342			break;
343	greg	2.6	}
344	greg	2.17	}
345	greg	2.6	}
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			}