src/util/rttree_reduce.c

#ifndef lint
static const char RCSid[] = "$Id: rttree_reduce.c,v 2.4 2011/06/03 18:12:58 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()
{
        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);
}

/* Load data array, filling zeroes for rank 3 demi-tensor */
static void
load_data()
{
        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");
}

/* Load BSDF array, coalesce uniform regions and format as tensor tree */
int
main(int argc, char *argv[])
{
        int     doheader = 1;
        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 '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], "rb", stdin) == NULL) {
                sprintf(errmsg, "cannot open input file '%s'", argv[i]);
                error(SYSTEM, errmsg);
        }
        if (infmt != 'a')
                SET_FILE_BINARY(stdin);
        load_data();
        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][-f{a|f|d}][-r {3|4}][-g log2grid][-t trim%%] [input]\n",
                        argv[0]);
        return(1);
}
Revision:	2.5
Committed:	Sat Aug 20 02:46:13 2011 UTC (12 years, 8 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R1
Changes since 2.4:	+4 -4 lines
Log Message:	Fixed bug in data ordering
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.5	static const char RCSid[] = "$Id: rttree_reduce.c,v 2.4 2011/06/03 18:12:58 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			set_threshold()
114			{
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			/* Load data array, filling zeroes for rank 3 demi-tensor */
239			static void
240			load_data()
241			{
242			int (readf)(float , int) = NULL;
243
244			switch (infmt) {
245			case 'a':
246			readf = &read_ascii;
247			break;
248			case 'f':
249			readf = &read_float;
250			break;
251			case 'd':
252			readf = &read_double;
253			break;
254			default:
255			error(COMMAND, "unsupported input format");
256			break;
257			}
258			datarr = (float )calloc(1<<(log2gttrank), sizeof(float));
259			if (datarr == NULL)
260			error(SYSTEM, "out of memory in load_data");
261			if (ttrank == 3) {
262			int ix, ox;
263	greg	2.4	for (ix = 0; ix < 1<<(log2g-1); ix++)
264	greg	2.1	for (ox = 0; ox < 1<<log2g; ox++)
265	greg	2.5	(*readf)(datarr+(((ix<<log2g)+ox)<<log2g),
266	greg	2.4	1<<log2g);
267	greg	2.1	} else /* ttrank == 4 */ {
268			int ix, iy, ox;
269			for (ix = 0; ix < 1<<log2g; ix++)
270			for (iy = 0; iy < 1<<log2g; iy++)
271			for (ox = 0; ox < 1<<log2g; ox++)
272			(*readf)(datarr +
273	greg	2.5	(((((ix<<log2g)+iy)<<log2g)+ox)<<log2g),
274	greg	2.1	1<<log2g);
275			}
276			(readf)(NULL, 0); / releases any buffers */
277			if (infmt == 'a') {
278			int c;
279			while ((c = getc(stdin)) != EOF) {
280			switch (c) {
281			case ' ':
282			case '\t':
283			case '\r':
284			case '\n':
285			continue;
286			}
287			error(WARNING, "data past end of expected input");
288			break;
289			}
290			} else if (getc(stdin) != EOF)
291			error(WARNING, "binary data past end of expected input");
292			}
293
294			/* Load BSDF array, coalesce uniform regions and format as tensor tree */
295			int
296			main(int argc, char *argv[])
297			{
298			int doheader = 1;
299			int bmin[4];
300			TNODE gtree;
301			int i;
302			/* get options and parameters */
303			for (i = 1; i < argc && argv[i][0] == '-'; i++)
304			switch (argv[i][1]) {
305			case 'h':
306			doheader = !doheader;
307			break;
308			case 'r':
309			ttrank = atoi(argv[++i]);
310			if (ttrank != 3 && ttrank != 4)
311			goto userr;
312			break;
313			case 'g':
314			log2g = atoi(argv[++i]);
315			if (log2g <= 1)
316			goto userr;
317			break;
318			case 't':
319	greg	2.2	pctcull = atof(argv[++i]);
320			if ((pctcull < 0) \| (pctcull >= 100.))
321	greg	2.1	goto userr;
322			break;
323			case 'f':
324			infmt = argv[i][2];
325			if (!infmt \|\| strchr("afd", infmt) == NULL)
326			goto userr;
327			break;
328			default:
329			goto userr;
330			}
331			if (i < argc-1)
332			goto userr;
333			/* load input data */
334			if (i == argc-1 && freopen(argv[i], "rb", stdin) == NULL) {
335			sprintf(errmsg, "cannot open input file '%s'", argv[i]);
336			error(SYSTEM, errmsg);
337			}
338			if (infmt != 'a')
339			SET_FILE_BINARY(stdin);
340			load_data();
341			if (doheader) {
342			for (i = 0; i < argc; i++) {
343			fputs(argv[i], stdout);
344			fputc(i < argc-1 ? ' ' : '\n', stdout);
345			}
346			fputc('\n', stdout);
347			}
348			gtree.kid = NULL; /* create our tree */
349			bmin[0] = bmin[1] = bmin[2] = bmin[3] = 0;
350			build_tree(&gtree, bmin, log2g);
351	greg	2.2	/* compute threshold & trim tree */
352			set_threshold();
353			trim_tree(&gtree);
354	greg	2.1	/* format to stdout */
355			print_tree(&gtree, bmin, log2g);
356			/* Clean up isn't necessary for main()...
357			free_kids(&gtree);
358			free(datarr);
359			*/
360			return(0);
361			userr:
362	greg	2.2	fprintf(stderr, "Usage: %s [-h][-f{a\|f\|d}][-r {3\|4}][-g log2grid][-t trim%%] [input]\n",
363	greg	2.1	argv[0]);
364			return(1);
365			}