ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/radiance/ray/src/common/bsdf_t.c

Comparing ray/src/common/bsdf_t.c (file contents):
Revision 3.1 by greg, Fri Feb 18 00:40:25 2011 UTC vs.
Revision 3.54 by greg, Tue Jan 25 01:34:20 2022 UTC

#	Line 1 | Line 1
1	+	#ifndef lint
2	+	static const char RCSid[] = "$Id$";
3	+	#endif
4		/*
5		*  bsdf_t.c
6		*
#	Line 7 | Line 10
10		*
11		*/
12
13	<	#include <stdio.h>
13	>	#define _USE_MATH_DEFINES
14	>	#include "rtio.h"
15		#include <stdlib.h>
16	+	#include <math.h>
17	+	#include <ctype.h>
18		#include "ezxml.h"
19		#include "bsdf.h"
20		#include "bsdf_t.h"
21	+	#include "hilbert.h"
22
23	+	/* Callback function type for SDtraverseTre() */
24	+	typedef int     SDtreCallback(float val, const double *cmin, double csiz,
25	+	void *cptr);
26	+	/* reference width maximum (1.0) */
27	+	static const unsigned   iwbits = sizeof(unsigned)*4;
28	+	static const unsigned   iwmax = 1<<(sizeof(unsigned)*4);
29	+	/* maximum cumulative value */
30	+	static const unsigned   cumlmax = ~0;
31	+	/* constant z-vector */
32	+	static const FVECT      zvec = {.0, .0, 1.};
33	+	/* quantization value */
34	+	static double           quantum = 1./256.;
35	+	/* our RGB primaries */
36	+	static C_COLOR          tt_RGB_prim[3];
37	+	static float            tt_RGB_coef[3];
38	+
39	+	static const double     czero[SD_MAXDIM];
40	+
41	+	enum {tt_Y, tt_u, tt_v};                /* tree components (tt_Y==0) */
42	+
43	+	/* Struct used for our distribution-building callback */
44	+	typedef struct {
45	+	short           nic;            /* number of input coordinates */
46	+	short           rev;            /* reversing query */
47	+	unsigned        alen;           /* current array length */
48	+	unsigned        nall;           /* number of allocated entries */
49	+	unsigned        wmin;           /* minimum square size so far */
50	+	unsigned        wmax;           /* maximum square size */
51	+	struct outdir_s {
52	+	unsigned        hent;           /* entering Hilbert index */
53	+	int             wid;            /* this square size */
54	+	float           bsdf;           /* BSDF for this square */
55	+	}               *darr;          /* output direction array */
56	+	} SDdistScaffold;
57	+
58		/* Allocate a new scattering distribution node */
59		static SDNode *
60		SDnewNode(int nd, int lg)
#	Line 30 | Line 72 | SDnewNode(int nd, int lg)
72		}
73		if (lg < 0) {
74		st = (SDNode *)malloc(sizeof(SDNode) +
75	<	((1<<nd) - 1)*sizeof(st->u.t[0]));
76	<	if (st != NULL)
77	<	memset(st->u.t, 0, (1<<nd)*sizeof(st->u.t[0]));
78	<	} else
75	>	sizeof(st->u.t[0])*((1<<nd) - 1));
76	>	if (st == NULL) {
77	>	sprintf(SDerrorDetail,
78	>	"Cannot allocate %d branch BSDF tree", 1<<nd);
79	>	return NULL;
80	>	}
81	>	memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
82	>	} else {
83		st = (SDNode *)malloc(sizeof(SDNode) +
84	<	((1 << nd*lg) - 1)*sizeof(st->u.v[0]));
85	<
40	<	if (st == NULL) {
41	<	if (lg < 0)
84	>	sizeof(st->u.v[0])*((1 << nd*lg) - 1));
85	>	if (st == NULL) {
86		sprintf(SDerrorDetail,
43	–	"Cannot allocate %d branch BSDF tree", nd);
44	–	else
45	–	sprintf(SDerrorDetail,
87		"Cannot allocate %d BSDF leaves", 1 << nd*lg);
88	<	return NULL;
88	>	return NULL;
89	>	}
90		}
91		st->ndim = nd;
92		st->log2GR = lg;
#	Line 53 | Line 95 | SDnewNode(int nd, int lg)
95
96		/* Free an SD tree */
97		static void
98	<	SDfreeTree(void *p)
98	>	SDfreeTre(SDNode *st)
99		{
100	<	SDNode  *st = (SDNode *)p;
59	<	int     i;
100	>	int     n;
101
102		if (st == NULL)
103		return;
104	<	for (i = (st->log2GR < 0) << st->ndim; i--; )
105	<	SDfreeTree(st->u.t[i]);
106	<	free((void *)st);
104	>	for (n = (st->log2GR < 0) << st->ndim; n--; )
105	>	SDfreeTre(st->u.t[n]);
106	>	free(st);
107		}
108
109	+	/* Free a variable-resolution BSDF */
110	+	static void
111	+	SDFreeBTre(void *p)
112	+	{
113	+	SDTre   *sdt = (SDTre *)p;
114	+
115	+	if (sdt == NULL)
116	+	return;
117	+	SDfreeTre(sdt->stc[tt_Y]);
118	+	SDfreeTre(sdt->stc[tt_u]);
119	+	SDfreeTre(sdt->stc[tt_v]);
120	+	free(sdt);
121	+	}
122	+
123	+	/* Fill branch's worth of grid values from subtree */
124	+	static void
125	+	fill_grid_branch(float *dptr, const float *sptr, int nd, int shft)
126	+	{
127	+	unsigned        n = 1 << (shft-1);
128	+
129	+	if (!--nd) {                    /* end on the line */
130	+	memcpy(dptr, sptr, sizeof(*dptr)*n);
131	+	return;
132	+	}
133	+	while (n--)                     /* recurse on each slice */
134	+	fill_grid_branch(dptr + (n << shft*nd),
135	+	sptr + (n << (shft-1)*nd), nd, shft);
136	+	}
137	+
138	+	/* Get pointer at appropriate offset for the given branch */
139	+	static float *
140	+	grid_branch_start(SDNode *st, int n)
141	+	{
142	+	unsigned        skipsiz = 1 << (st->log2GR - 1);
143	+	float           *vptr = st->u.v;
144	+	int             i;
145	+
146	+	for (i = st->ndim; i--; skipsiz <<= st->log2GR)
147	+	if (1<<i & n)
148	+	vptr += skipsiz;
149	+	return vptr;
150	+	}
151	+
152	+	/* Simplify (consolidate) a tree by flattening uniform depth regions */
153	+	static SDNode *
154	+	SDsimplifyTre(SDNode *st)
155	+	{
156	+	int             match, n;
157	+
158	+	if (st == NULL)                 /* check for invalid tree */
159	+	return NULL;
160	+	if (st->log2GR >= 0)            /* grid just returns unaltered */
161	+	return st;
162	+	match = 1;                      /* check if grids below match */
163	+	for (n = 0; n < 1<<st->ndim; n++) {
164	+	if ((st->u.t[n] = SDsimplifyTre(st->u.t[n])) == NULL)
165	+	return NULL;    /* propogate error up call stack */
166	+	match &= (st->u.t[n]->log2GR == st->u.t[0]->log2GR);
167	+	}
168	+	if (match && (match = st->u.t[0]->log2GR) >= 0) {
169	+	SDNode  *stn = SDnewNode(st->ndim, match + 1);
170	+	if (stn == NULL)        /* out of memory? */
171	+	return st;
172	+	/* transfer values to new grid */
173	+	for (n = 1 << st->ndim; n--; )
174	+	fill_grid_branch(grid_branch_start(stn, n),
175	+	st->u.t[n]->u.v, stn->ndim, stn->log2GR);
176	+	SDfreeTre(st);          /* free old tree */
177	+	st = stn;               /* return new one */
178	+	}
179	+	return st;
180	+	}
181	+
182	+	/* Assign the given voxel in tree (produces no grid nodes) */
183	+	static SDNode *
184	+	SDsetVoxel(SDNode *sroot, int nd, const double *tmin, const double tsiz, float val)
185	+	{
186	+	double  ctrk[SD_MAXDIM];
187	+	double  csiz = 1.;
188	+	SDNode  *st;
189	+	int     i, n;
190	+	/* check arguments */
191	+	for (i = nd; i-- > 0; )
192	+	if ((tmin[i] < .0) | (tmin[i] >= 1.-FTINY))
193	+	break;
194	+	if ((i >= 0) | (nd <= 0) | (tsiz <= FTINY) | (tsiz > 1.+FTINY) |
195	+	(sroot != NULL && sroot->ndim != nd)) {
196	+	SDfreeTre(sroot);
197	+	return NULL;
198	+	}
199	+	if (tsiz >= 1.-FTINY) {         /* special case when tree is a leaf */
200	+	SDfreeTre(sroot);
201	+	if ((sroot = SDnewNode(nd, 0)) != NULL)
202	+	sroot->u.v[0] = val;
203	+	return sroot;
204	+	}
205	+	/* make sure we have branching root */
206	+	if (sroot != NULL && sroot->log2GR >= 0) {
207	+	SDfreeTre(sroot); sroot = NULL;
208	+	}
209	+	if (sroot == NULL && (sroot = SDnewNode(nd, -1)) == NULL)
210	+	return NULL;
211	+	st = sroot;                     /* climb/grow tree */
212	+	memset(ctrk, 0, sizeof(ctrk));
213	+	for ( ; ; ) {
214	+	csiz *= .5;             /* find appropriate branch */
215	+	n = 0;
216	+	for (i = nd; i--; )
217	+	if (ctrk[i]+csiz <= tmin[i]+FTINY) {
218	+	ctrk[i] += csiz;
219	+	n |= 1 << i;
220	+	}
221	+	/* reached desired voxel? */
222	+	if (csiz <= tsiz+FTINY) {
223	+	SDfreeTre(st->u.t[n]);
224	+	st = st->u.t[n] = SDnewNode(nd, 0);
225	+	break;
226	+	}
227	+	/* else grow tree as needed */
228	+	if (st->u.t[n] != NULL && st->u.t[n]->log2GR >= 0) {
229	+	SDfreeTre(st->u.t[n]); st->u.t[n] = NULL;
230	+	}
231	+	if (st->u.t[n] == NULL)
232	+	st->u.t[n] = SDnewNode(nd, -1);
233	+	if ((st = st->u.t[n]) == NULL)
234	+	break;
235	+	}
236	+	if (st == NULL) {
237	+	SDfreeTre(sroot);
238	+	return NULL;
239	+	}
240	+	st->u.v[0] = val;               /* assign leaf and return root */
241	+	return sroot;
242	+	}
243	+
244	+	/* Find smallest leaf in tree */
245	+	static double
246	+	SDsmallestLeaf(const SDNode *st)
247	+	{
248	+	if (st->log2GR < 0) {           /* tree branches */
249	+	double  lmin = 1.;
250	+	int     n;
251	+	for (n = 1<<st->ndim; n--; ) {
252	+	double  lsiz = SDsmallestLeaf(st->u.t[n]);
253	+	if (lsiz < lmin)
254	+	lmin = lsiz;
255	+	}
256	+	return .5*lmin;
257	+	}
258	+	/* leaf grid width */
259	+	return 1. / (double)(1 << st->log2GR);
260	+	}
261	+
262		/* Add up N-dimensional hypercube array values over the given box */
263		static double
264	<	SDiterSum(const float *va, int nd, int siz, const int *imin, const int *imax)
264	>	SDiterSum(const float *va, int nd, int shft, const int *imin, const int *imax)
265		{
266	+	const unsigned  skipsiz = 1 << --nd*shft;
267		double          sum = .0;
73	–	unsigned        skipsiz = 1;
268		int             i;
269	<
270	<	for (i = nd; --i > 0; )
271	<	skipsiz *= siz;
269	>
270	>	va += *imin * skipsiz;
271	>
272		if (skipsiz == 1)
273		for (i = *imin; i < *imax; i++)
274	<	sum += va[i];
274	>	sum += *va++;
275		else
276	<	for (i = *imin; i < *imax; i++)
277	<	sum += SDiterSum(va + i*skipsiz,
84	<	nd-1, siz, imin+1, imax+1);
276	>	for (i = *imin; i < *imax; i++, va += skipsiz)
277	>	sum += SDiterSum(va, nd, shft, imin+1, imax+1);
278		return sum;
279		}
280
281		/* Average BSDF leaves over an orthotope defined by the unit hypercube */
282		static double
283	<	SDavgBox(const SDNode *st, const double *bmin, const double *bmax)
283	>	SDavgTreBox(const SDNode *st, const double *bmin, const double *bmax)
284		{
92	–	int             imin[SD_MAXDIM], imax[SD_MAXDIM];
285		unsigned        n;
286		int             i;
287
#	Line 99 | Line 291 | SDavgBox(const SDNode *st, const double *bmin, const d
291		for (i = st->ndim; i--; ) {
292		if (bmin[i] >= 1.)
293		return .0;
294	<	if (bmax[i] <= .0)
294	>	if (bmax[i] <= 0)
295		return .0;
296		if (bmin[i] >= bmax[i])
297		return .0;
#	Line 107 | Line 299 | SDavgBox(const SDNode *st, const double *bmin, const d
299		if (st->log2GR < 0) {           /* iterate on subtree */
300		double          sum = .0, wsum = 1e-20;
301		double          sbmin[SD_MAXDIM], sbmax[SD_MAXDIM], w;
110	–
302		for (n = 1 << st->ndim; n--; ) {
303		w = 1.;
304		for (i = st->ndim; i--; ) {
#	Line 119 | Line 310 | SDavgBox(const SDNode *st, const double *bmin, const d
310		}
311		if (sbmin[i] < .0) sbmin[i] = .0;
312		if (sbmax[i] > 1.) sbmax[i] = 1.;
313	+	if (sbmin[i] >= sbmax[i]) {
314	+	w = .0;
315	+	break;
316	+	}
317		w *= sbmax[i] - sbmin[i];
318		}
319		if (w > 1e-10) {
320	<	sum += w * SDavgBox(st->u.t[n], sbmin, sbmax);
320	>	sum += w * SDavgTreBox(st->u.t[n], sbmin, sbmax);
321		wsum += w;
322		}
323		}
324		return sum / wsum;
325	+	} else {                        /* iterate over leaves */
326	+	int             imin[SD_MAXDIM], imax[SD_MAXDIM];
327	+
328	+	n = 1;
329	+	for (i = st->ndim; i--; ) {
330	+	imin[i] = (bmin[i] <= 0) ? 0 :
331	+	(int)((1 << st->log2GR)*bmin[i]);
332	+	imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR) :
333	+	(int)((1 << st->log2GR)*bmax[i] + .999999);
334	+	n *= imax[i] - imin[i];
335	+	}
336	+	if (n)
337	+	return SDiterSum(st->u.v, st->ndim,
338	+	st->log2GR, imin, imax) / (double)n;
339		}
340	<	n = 1;                          /* iterate over leaves */
341	<	for (i = st->ndim; i--; ) {
342	<	imin[i] = (bmin[i] <= .0) ? 0
343	<	: (int)((1 << st->log2GR)*bmin[i]);
344	<	imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR)
345	<	: (int)((1 << st->log2GR)*bmax[i] + .999999);
346	<	n *= imax[i] - imin[i];
340	>	return .0;
341	>	}
342	>
343	>	/* Recursive call for SDtraverseTre() */
344	>	static int
345	>	SDdotravTre(const SDNode *st, const double *pos, int cmask,
346	>	SDtreCallback *cf, void *cptr,
347	>	const double *cmin, double csiz)
348	>	{
349	>	int     rv, rval = 0;
350	>	double  bmin[SD_MAXDIM];
351	>	int     i, n;
352	>	/* paranoia */
353	>	if (st == NULL)
354	>	return 0;
355	>	/* in branches? */
356	>	if (st->log2GR < 0) {
357	>	unsigned        skipmask = 0;
358	>	csiz *= .5;
359	>	for (i = st->ndim; i--; )
360	>	if (1<<i & cmask) {
361	>	if (pos[i] < cmin[i] + csiz)
362	>	for (n = 1 << st->ndim; n--; ) {
363	>	if (n & 1<<i)
364	>	skipmask |= 1<<n;
365	>	}
366	>	else
367	>	for (n = 1 << st->ndim; n--; ) {
368	>	if (!(n & 1<<i))
369	>	skipmask |= 1<<n;
370	>	}
371	>	}
372	>	for (n = 1 << st->ndim; n--; ) {
373	>	if (1<<n & skipmask)
374	>	continue;
375	>	for (i = st->ndim; i--; )
376	>	if (1<<i & n)
377	>	bmin[i] = cmin[i] + csiz;
378	>	else
379	>	bmin[i] = cmin[i];
380	>
381	>	rval += rv = SDdotravTre(st->u.t[n], pos, cmask,
382	>	cf, cptr, bmin, csiz);
383	>	if (rv < 0)
384	>	return rv;
385	>	}
386	>	} else {                        /* else traverse leaves */
387	>	int     clim[SD_MAXDIM][2];
388	>	int     cpos[SD_MAXDIM];
389	>
390	>	if (st->log2GR == 0)    /* short cut */
391	>	return (*cf)(st->u.v[0], cmin, csiz, cptr);
392	>
393	>	csiz /= (double)(1 << st->log2GR);
394	>	/* assign coord. ranges */
395	>	for (i = st->ndim; i--; )
396	>	if (1<<i & cmask) {
397	>	clim[i][0] = (pos[i] - cmin[i])/csiz;
398	>	/* check overflow from f.p. error */
399	>	clim[i][0] -= clim[i][0] >> st->log2GR;
400	>	clim[i][1] = clim[i][0] + 1;
401	>	} else {
402	>	clim[i][0] = 0;
403	>	clim[i][1] = 1 << st->log2GR;
404	>	}
405	>	#if (SD_MAXDIM == 4)
406	>	bmin[0] = cmin[0] + csiz*clim[0][0];
407	>	for (cpos[0] = clim[0][0]; cpos[0] < clim[0][1]; cpos[0]++) {
408	>	bmin[1] = cmin[1] + csiz*clim[1][0];
409	>	for (cpos[1] = clim[1][0]; cpos[1] < clim[1][1]; cpos[1]++) {
410	>	bmin[2] = cmin[2] + csiz*clim[2][0];
411	>	if (st->ndim == 3) {
412	>	cpos[2] = clim[2][0];
413	>	n = cpos[0];
414	>	for (i = 1; i < 3; i++)
415	>	n = (n << st->log2GR) + cpos[i];
416	>	for ( ; cpos[2] < clim[2][1]; cpos[2]++) {
417	>	rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr);
418	>	if (rv < 0)
419	>	return rv;
420	>	bmin[2] += csiz;
421	>	}
422	>	} else {
423	>	for (cpos[2] = clim[2][0]; cpos[2] < clim[2][1]; cpos[2]++) {
424	>	bmin[3] = cmin[3] + csiz*(cpos[3] = clim[3][0]);
425	>	n = cpos[0];
426	>	for (i = 1; i < 4; i++)
427	>	n = (n << st->log2GR) + cpos[i];
428	>	for ( ; cpos[3] < clim[3][1]; cpos[3]++) {
429	>	rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr);
430	>	if (rv < 0)
431	>	return rv;
432	>	bmin[3] += csiz;
433	>	}
434	>	bmin[2] += csiz;
435	>	}
436	>	}
437	>	bmin[1] += csiz;
438	>	}
439	>	bmin[0] += csiz;
440	>	}
441	>	#else
442	>	_!_ "broken code segment!"
443	>	#endif
444		}
445	<	if (!n)
445	>	return rval;
446	>	}
447	>
448	>	/* Traverse a tree, visiting nodes in a slice that fits partial position */
449	>	static int
450	>	SDtraverseTre(const SDNode *st, const double *pos, int cmask,
451	>	SDtreCallback *cf, void *cptr)
452	>	{
453	>	int             i;
454	>	/* check arguments */
455	>	if ((st == NULL) | (cf == NULL))
456	>	return -1;
457	>	for (i = st->ndim; i--; )
458	>	if (1<<i & cmask && (pos[i] < 0) | (pos[i] >= 1.))
459	>	return -1;
460	>
461	>	return SDdotravTre(st, pos, cmask, cf, cptr, czero, 1.);
462	>	}
463	>
464	>	/* Look up tree value at the given grid position */
465	>	static float
466	>	SDlookupTre(const SDNode *st, const double *pos, double *hcube)
467	>	{
468	>	double  spos[SD_MAXDIM];
469	>	int     i, n, t;
470	>	/* initialize voxel return */
471	>	if (hcube) {
472	>	hcube[i = st->ndim] = 1.;
473	>	while (i--)
474	>	hcube[i] = .0;
475	>	}
476	>	/* climb the tree */
477	>	while (st != NULL && st->log2GR < 0) {
478	>	n = 0;                  /* move to appropriate branch */
479	>	if (hcube) hcube[st->ndim] *= .5;
480	>	for (i = st->ndim; i--; ) {
481	>	spos[i] = 2.*pos[i];
482	>	t = (spos[i] >= 1.);
483	>	n |= t<<i;
484	>	spos[i] -= (double)t;
485	>	if (hcube) hcube[i] += (double)t * hcube[st->ndim];
486	>	}
487	>	st = st->u.t[n];        /* avoids tail recursion */
488	>	pos = spos;
489	>	}
490	>	if (st == NULL)                 /* should never happen? */
491		return .0;
492	<
493	<	return SDiterSum(st->u.v, st->ndim, 1 << st->log2GR, imin, imax) /
494	<	(double)n;
492	>	if (st->log2GR == 0)            /* short cut */
493	>	return st->u.v[0];
494	>	n = t = 0;                      /* find grid array index */
495	>	for (i = st->ndim; i--; ) {
496	>	n += (int)((1<<st->log2GR)*pos[i]) << t;
497	>	t += st->log2GR;
498	>	}
499	>	if (hcube) {                    /* compute final hypercube */
500	>	hcube[st->ndim] /= (double)(1<<st->log2GR);
501	>	for (i = st->ndim; i--; )
502	>	hcube[i] += floor((1<<st->log2GR)*pos[i])*hcube[st->ndim];
503	>	}
504	>	return st->u.v[n];              /* no interpolation */
505		}
506
507	+	/* Convert CIE (Y,u',v') color to our RGB */
508	+	static void
509	+	SDyuv2rgb(double yval, double uprime, double vprime, float rgb[3])
510	+	{
511	+	const double    dfact = 1./(6.*uprime - 16.*vprime + 12.);
512	+	C_COLOR         cxy;
513	+
514	+	c_cset(&cxy, 9.*uprime*dfact, 4.*vprime*dfact);
515	+	c_toSharpRGB(&cxy, yval, rgb);
516	+	}
517	+
518	+	static double
519	+	pfrac(double x)
520	+	{
521	+	return( x - (int)x );
522	+	}
523	+
524	+	/* Query BSDF value and sample hypercube for the given vectors */
525	+	static int
526	+	SDqueryTre(const SDTre *sdt, float *coef,
527	+	const FVECT inVec, const FVECT outVec, double *hc)
528	+	{
529	+	const RREAL     *vtmp;
530	+	double          hcube[SD_MAXDIM+1];
531	+	float           yval;
532	+	FVECT           rOutVec;
533	+	RREAL           gridPos[4];
534	+	double          d;
535	+	int             i;
536	+
537	+	if (sdt->stc[tt_Y] == NULL)     /* paranoia, I hope */
538	+	return 0;
539	+
540	+	switch (sdt->sidef) {           /* whose side are you on? */
541	+	case SD_FREFL:
542	+	if ((outVec[2] < 0) | (inVec[2] < 0))
543	+	return 0;
544	+	break;
545	+	case SD_BREFL:
546	+	if ((outVec[2] > 0) | (inVec[2] > 0))
547	+	return 0;
548	+	break;
549	+	case SD_FXMIT:
550	+	if (outVec[2] > 0) {
551	+	if (inVec[2] > 0)
552	+	return 0;
553	+	vtmp = outVec; outVec = inVec; inVec = vtmp;
554	+	} else if (inVec[2] < 0)
555	+	return 0;
556	+	break;
557	+	case SD_BXMIT:
558	+	if (inVec[2] > 0) {
559	+	if (outVec[2] > 0)
560	+	return 0;
561	+	vtmp = outVec; outVec = inVec; inVec = vtmp;
562	+	} else if (outVec[2] < 0)
563	+	return 0;
564	+	break;
565	+	default:
566	+	return 0;
567	+	}
568	+	/* convert vector coordinates */
569	+	if (sdt->stc[tt_Y]->ndim == 3) {
570	+	spinvector(rOutVec, outVec, zvec, -atan2(-inVec[1],-inVec[0]));
571	+	gridPos[0] = (.5-FTINY) -
572	+	.5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
573	+	disk2square(gridPos+1, rOutVec[0], rOutVec[1]);
574	+	} else if (sdt->stc[tt_Y]->ndim == 4) {
575	+	disk2square(gridPos, -inVec[0], -inVec[1]);
576	+	disk2square(gridPos+2, outVec[0], outVec[1]);
577	+	} else
578	+	return 0;               /* should be internal error */
579	+
580	+	if (hc == NULL) hc = hcube;     /* get BSDF value */
581	+	yval = SDlookupTre(sdt->stc[tt_Y], gridPos, hc);
582	+	if (coef == NULL)               /* just getting hypercube? */
583	+	return 1;
584	+	d = 0;                          /* position-specific perturbation */
585	+	for (i = sdt->stc[tt_Y]->ndim; i--; )
586	+	d += pfrac((2<<i)/(hc[i]+.01687)) - .5;
587	+	yval *= 1. + 1e-4*d;            /* assumes tolerance is > 0.04% */
588	+	if (sdt->stc[tt_u] == NULL || sdt->stc[tt_v] == NULL) {
589	+	*coef = yval;
590	+	return 1;               /* no color */
591	+	}
592	+	/* else decode color */
593	+	SDyuv2rgb(yval, SDlookupTre(sdt->stc[tt_u], gridPos, NULL),
594	+	SDlookupTre(sdt->stc[tt_v], gridPos, NULL), coef);
595	+	coef[0] *= tt_RGB_coef[0];
596	+	coef[1] *= tt_RGB_coef[1];
597	+	coef[2] *= tt_RGB_coef[2];
598	+	return 3;
599	+	}
600	+
601	+	/* Compute non-diffuse component for variable-resolution BSDF */
602	+	static int
603	+	SDgetTreBSDF(float coef[SDmaxCh], const FVECT outVec,
604	+	const FVECT inVec, SDComponent *sdc)
605	+	{
606	+	/* check arguments */
607	+	if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
608	+	|| sdc->dist == NULL)
609	+	return 0;
610	+	/* get nearest BSDF value */
611	+	return SDqueryTre((SDTre *)sdc->dist, coef, outVec, inVec, NULL);
612	+	}
613	+
614	+	/* Callback to build cumulative distribution using SDtraverseTre() */
615	+	static int
616	+	build_scaffold(float val, const double *cmin, double csiz, void *cptr)
617	+	{
618	+	SDdistScaffold  *sp = (SDdistScaffold *)cptr;
619	+	int             wid = csiz*(double)iwmax + .5;
620	+	double          revcmin[2];
621	+	bitmask_t       bmin[2], bmax[2];
622	+
623	+	if (sp->rev) {                  /* need to reverse sense? */
624	+	revcmin[0] = 1. - cmin[0] - csiz;
625	+	revcmin[1] = 1. - cmin[1] - csiz;
626	+	cmin = revcmin;
627	+	} else {
628	+	cmin += sp->nic;        /* else skip to output coords */
629	+	}
630	+	if (wid < sp->wmin)             /* new minimum width? */
631	+	sp->wmin = wid;
632	+	if (wid > sp->wmax)             /* new maximum? */
633	+	sp->wmax = wid;
634	+	if (sp->alen >= sp->nall) {     /* need more space? */
635	+	struct outdir_s *ndarr;
636	+	sp->nall = (int)(1.5*sp->nall) + 256;
637	+	ndarr = (struct outdir_s *)realloc(sp->darr,
638	+	sizeof(struct outdir_s)*sp->nall);
639	+	if (ndarr == NULL) {
640	+	sprintf(SDerrorDetail,
641	+	"Cannot grow scaffold to %u entries", sp->nall);
642	+	return -1;      /* abort build */
643	+	}
644	+	sp->darr = ndarr;
645	+	}
646	+	/* find Hilbert entry index */
647	+	bmin[0] = cmin[0]*(double)iwmax + .5;
648	+	bmin[1] = cmin[1]*(double)iwmax + .5;
649	+	bmax[0] = bmin[0] + wid-1;
650	+	bmax[1] = bmin[1] + wid-1;
651	+	hilbert_box_vtx(2, sizeof(bitmask_t), iwbits, 1, bmin, bmax);
652	+	sp->darr[sp->alen].hent = hilbert_c2i(2, iwbits, bmin);
653	+	sp->darr[sp->alen].wid = wid;
654	+	sp->darr[sp->alen].bsdf = val;
655	+	sp->alen++;                     /* on to the next entry */
656	+	return 0;
657	+	}
658	+
659	+	/* Scaffold comparison function for qsort -- ascending Hilbert index */
660	+	static int
661	+	sscmp(const void *p1, const void *p2)
662	+	{
663	+	unsigned        h1 = (*(const struct outdir_s *)p1).hent;
664	+	unsigned        h2 = (*(const struct outdir_s *)p2).hent;
665	+
666	+	if (h1 > h2)
667	+	return 1;
668	+	if (h1 < h2)
669	+	return -1;
670	+	return 0;
671	+	}
672	+
673	+	/* Create a new cumulative distribution for the given input direction */
674	+	static SDTreCDst *
675	+	make_cdist(const SDTre *sdt, const double *invec, int rev)
676	+	{
677	+	SDdistScaffold  myScaffold;
678	+	double          pos[4];
679	+	int             cmask;
680	+	SDTreCDst       *cd;
681	+	struct outdir_s *sp;
682	+	double          scale, cursum;
683	+	int             i;
684	+	/* initialize scaffold */
685	+	myScaffold.wmin = iwmax;
686	+	myScaffold.wmax = 0;
687	+	myScaffold.nic = sdt->stc[tt_Y]->ndim - 2;
688	+	myScaffold.rev = rev;
689	+	myScaffold.alen = 0;
690	+	myScaffold.nall = 512;
691	+	myScaffold.darr = (struct outdir_s *)malloc(sizeof(struct outdir_s) *
692	+	myScaffold.nall);
693	+	if (myScaffold.darr == NULL)
694	+	return NULL;
695	+	/* set up traversal */
696	+	cmask = (1<<myScaffold.nic) - 1;
697	+	for (i = myScaffold.nic; i--; )
698	+	pos[i+2*rev] = invec[i];
699	+	cmask <<= 2*rev;
700	+	/* grow the distribution */
701	+	if (SDtraverseTre(sdt->stc[tt_Y], pos, cmask,
702	+	build_scaffold, &myScaffold) < 0) {
703	+	free(myScaffold.darr);
704	+	return NULL;
705	+	}
706	+	/* allocate result holder */
707	+	cd = (SDTreCDst *)malloc(sizeof(SDTreCDst) +
708	+	sizeof(cd->carr[0])*myScaffold.alen);
709	+	if (cd == NULL) {
710	+	sprintf(SDerrorDetail,
711	+	"Cannot allocate %u entry cumulative distribution",
712	+	myScaffold.alen);
713	+	free(myScaffold.darr);
714	+	return NULL;
715	+	}
716	+	cd->isodist = (myScaffold.nic == 1);
717	+	/* sort the distribution */
718	+	qsort(myScaffold.darr, cd->calen = myScaffold.alen,
719	+	sizeof(struct outdir_s), sscmp);
720	+
721	+	/* record input range */
722	+	scale = myScaffold.wmin / (double)iwmax;
723	+	for (i = myScaffold.nic; i--; ) {
724	+	cd->clim[i][0] = floor(pos[i+2*rev]/scale) * scale;
725	+	cd->clim[i][1] = cd->clim[i][0] + scale;
726	+	}
727	+	if (cd->isodist) {              /* avoid issue in SDqueryTreProjSA() */
728	+	cd->clim[1][0] = cd->clim[0][0];
729	+	cd->clim[1][1] = cd->clim[0][1];
730	+	}
731	+	cd->max_psa = myScaffold.wmax / (double)iwmax;
732	+	cd->max_psa *= cd->max_psa * M_PI;
733	+	if (rev)
734	+	cd->sidef = (sdt->sidef==SD_BXMIT) ? SD_FXMIT : SD_BXMIT;
735	+	else
736	+	cd->sidef = sdt->sidef;
737	+	cd->cTotal = 1e-20;             /* compute directional total */
738	+	sp = myScaffold.darr;
739	+	for (i = myScaffold.alen; i--; sp++)
740	+	cd->cTotal += sp->bsdf * (double)sp->wid * sp->wid;
741	+	cursum = .0;                    /* go back and get cumulative values */
742	+	scale = (double)cumlmax / cd->cTotal;
743	+	sp = myScaffold.darr;
744	+	for (i = 0; i < cd->calen; i++, sp++) {
745	+	cd->carr[i].hndx = sp->hent;
746	+	cd->carr[i].cuml = scale*cursum + .5;
747	+	cursum += sp->bsdf * (double)sp->wid * sp->wid;
748	+	}
749	+	cd->carr[i].hndx = ~0;          /* make final entry */
750	+	cd->carr[i].cuml = cumlmax;
751	+	cd->cTotal *= M_PI/(double)iwmax/iwmax;
752	+	/* all done, clean up and return */
753	+	free(myScaffold.darr);
754	+	return cd;
755	+	}
756	+
757	+	/* Find or allocate a cumulative distribution for the given incoming vector */
758	+	const SDCDst *
759	+	SDgetTreCDist(const FVECT inVec, SDComponent *sdc)
760	+	{
761	+	unsigned long   cacheLeft = SDmaxCache;
762	+	const SDTre     *sdt;
763	+	RREAL           inCoord[2];
764	+	int             i;
765	+	int             mode;
766	+	SDTreCDst       *cd, *cdlast, *cdlimit;
767	+	/* check arguments */
768	+	if ((inVec == NULL) | (sdc == NULL) ||
769	+	(sdt = (SDTre *)sdc->dist) == NULL)
770	+	return NULL;
771	+	switch (mode = sdt->sidef) {    /* check direction */
772	+	case SD_FREFL:
773	+	if (inVec[2] < 0)
774	+	return NULL;
775	+	break;
776	+	case SD_BREFL:
777	+	if (inVec[2] > 0)
778	+	return NULL;
779	+	break;
780	+	case SD_FXMIT:
781	+	if (inVec[2] < 0)
782	+	mode = SD_BXMIT;
783	+	break;
784	+	case SD_BXMIT:
785	+	if (inVec[2] > 0)
786	+	mode = SD_FXMIT;
787	+	break;
788	+	default:
789	+	return NULL;
790	+	}
791	+	if (sdt->stc[tt_Y]->ndim == 3) {        /* isotropic BSDF? */
792	+	if (mode != sdt->sidef) /* XXX unhandled reciprocity */
793	+	return &SDemptyCD;
794	+	inCoord[0] = (.5-FTINY) -
795	+	.5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
796	+	} else if (sdt->stc[tt_Y]->ndim == 4) {
797	+	if (mode != sdt->sidef) /* use reciprocity? */
798	+	disk2square(inCoord, inVec[0], inVec[1]);
799	+	else
800	+	disk2square(inCoord, -inVec[0], -inVec[1]);
801	+	} else
802	+	return NULL;            /* should be internal error */
803	+	/* quantize to avoid f.p. errors */
804	+	for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
805	+	inCoord[i] = floor(inCoord[i]/quantum)*quantum + .5*quantum;
806	+	cdlast = cdlimit = NULL;        /* check for direction in cache list */
807	+	/* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
808	+	for (cd = (SDTreCDst *)sdc->cdList; cd != NULL;
809	+	cdlast = cd, cd = cd->next) {
810	+	if (cacheLeft) {        /* check cache size limit */
811	+	long    csiz = sizeof(SDTreCDst) +
812	+	sizeof(cd->carr[0])*cd->calen;
813	+	if (cacheLeft > csiz)
814	+	cacheLeft -= csiz;
815	+	else {
816	+	cdlimit = cdlast;
817	+	cacheLeft = 0;
818	+	}
819	+	}
820	+	if (cd->sidef != mode)
821	+	continue;
822	+	for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
823	+	if ((cd->clim[i][0] > inCoord[i]) |
824	+	(inCoord[i] >= cd->clim[i][1]))
825	+	break;
826	+	if (i < 0)
827	+	break;          /* means we have a match */
828	+	}
829	+	if (cd == NULL) {               /* need to create new entry? */
830	+	if (cdlimit != NULL)    /* exceeded cache size limit? */
831	+	while ((cd = cdlimit->next) != NULL) {
832	+	cdlimit->next = cd->next;
833	+	free(cd);
834	+	}
835	+	cdlast = cd = make_cdist(sdt, inCoord, mode != sdt->sidef);
836	+	}
837	+	if (cdlast != NULL) {           /* move entry to head of cache list */
838	+	cdlast->next = cd->next;
839	+	cd->next = (SDTreCDst *)sdc->cdList;
840	+	sdc->cdList = (SDCDst *)cd;
841	+	}
842	+	/* END MUTEX LOCK */
843	+	return (SDCDst *)cd;            /* ready to go */
844	+	}
845	+
846	+	/* Query solid angle for vector(s) */
847	+	static SDError
848	+	SDqueryTreProjSA(double *psa, const FVECT v1, const RREAL *v2,
849	+	int qflags, SDComponent *sdc)
850	+	{
851	+	double          myPSA[2];
852	+	/* check arguments */
853	+	if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
854	+	sdc->dist == NULL)
855	+	return SDEargument;
856	+	/* get projected solid angle(s) */
857	+	if (v2 != NULL) {
858	+	const SDTre     *sdt = (SDTre *)sdc->dist;
859	+	double          hcube[SD_MAXDIM+1];
860	+	if (!SDqueryTre(sdt, NULL, v1, v2, hcube)) {
861	+	strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA");
862	+	return SDEinternal;
863	+	}
864	+	myPSA[0] = hcube[sdt->stc[tt_Y]->ndim];
865	+	myPSA[1] = myPSA[0] *= myPSA[0] * M_PI;
866	+	} else {
867	+	const SDTreCDst *cd = (const SDTreCDst *)SDgetTreCDist(v1, sdc);
868	+	if (cd == NULL)
869	+	myPSA[0] = myPSA[1] = 0;
870	+	else {
871	+	myPSA[0] = M_PI * (cd->clim[0][1] - cd->clim[0][0]) *
872	+	(cd->clim[1][1] - cd->clim[1][0]);
873	+	myPSA[1] = cd->max_psa;
874	+	}
875	+	}
876	+	switch (qflags) {               /* record based on flag settings */
877	+	case SDqueryVal:
878	+	*psa = myPSA[0];
879	+	break;
880	+	case SDqueryMax:
881	+	if (myPSA[1] > *psa)
882	+	*psa = myPSA[1];
883	+	break;
884	+	case SDqueryMin+SDqueryMax:
885	+	if (myPSA[1] > psa[1])
886	+	psa[1] = myPSA[1];
887	+	/* fall through */
888	+	case SDqueryMin:
889	+	if ((myPSA[0] > 0) & (myPSA[0] < psa[0]))
890	+	psa[0] = myPSA[0];
891	+	break;
892	+	}
893	+	return SDEnone;
894	+	}
895	+
896	+	/* Sample cumulative distribution */
897	+	static SDError
898	+	SDsampTreCDist(FVECT ioVec, double randX, const SDCDst *cdp)
899	+	{
900	+	const unsigned  nBitsC = 4*sizeof(bitmask_t);
901	+	const unsigned  nExtraBits = 8*(sizeof(bitmask_t)-sizeof(unsigned));
902	+	const SDTreCDst *cd = (const SDTreCDst *)cdp;
903	+	const unsigned  target = randX*cumlmax;
904	+	bitmask_t       hndx, hcoord[2];
905	+	FVECT           gpos;
906	+	double          rotangle;
907	+	int             i, iupper, ilower;
908	+	/* check arguments */
909	+	if ((ioVec == NULL) | (cd == NULL))
910	+	return SDEargument;
911	+	if (!cd->sidef)
912	+	return SDEnone;         /* XXX should never happen */
913	+	if (ioVec[2] > 0) {
914	+	if ((cd->sidef != SD_FREFL) & (cd->sidef != SD_FXMIT))
915	+	return SDEargument;
916	+	} else if ((cd->sidef != SD_BREFL) & (cd->sidef != SD_BXMIT))
917	+	return SDEargument;
918	+	/* binary search to find position */
919	+	ilower = 0; iupper = cd->calen;
920	+	while ((i = (iupper + ilower) >> 1) != ilower)
921	+	if (target >= cd->carr[i].cuml)
922	+	ilower = i;
923	+	else
924	+	iupper = i;
925	+	/* localize random position */
926	+	randX = (randX*cumlmax - cd->carr[ilower].cuml) /
927	+	(double)(cd->carr[iupper].cuml - cd->carr[ilower].cuml);
928	+	/* index in longer Hilbert curve */
929	+	hndx = (randX*cd->carr[iupper].hndx + (1.-randX)*cd->carr[ilower].hndx)
930	+	* (double)((bitmask_t)1 << nExtraBits);
931	+	/* convert Hilbert index to vector */
932	+	hilbert_i2c(2, nBitsC, hndx, hcoord);
933	+	for (i = 2; i--; )
934	+	gpos[i] = ((double)hcoord[i] + rand()*(1./(RAND_MAX+.5))) /
935	+	(double)((bitmask_t)1 << nBitsC);
936	+	square2disk(gpos, gpos[0], gpos[1]);
937	+	/* compute Z-coordinate */
938	+	gpos[2] = 1. - gpos[0]*gpos[0] - gpos[1]*gpos[1];
939	+	gpos[2] = sqrt(gpos[2]*(gpos[2]>0));
940	+	/* emit from back? */
941	+	if ((cd->sidef == SD_BREFL) | (cd->sidef == SD_FXMIT))
942	+	gpos[2] = -gpos[2];
943	+	if (cd->isodist) {              /* rotate isotropic sample */
944	+	rotangle = atan2(-ioVec[1],-ioVec[0]);
945	+	spinvector(ioVec, gpos, zvec, rotangle);
946	+	} else
947	+	VCOPY(ioVec, gpos);
948	+	return SDEnone;
949	+	}
950	+
951	+	/* Advance pointer to the next non-white character in the string (or nul) */
952	+	static int
953	+	next_token(char **spp)
954	+	{
955	+	while (isspace(**spp))
956	+	++*spp;
957	+	return **spp;
958	+	}
959	+
960	+	/* Advance pointer past matching token (or any token if c==0) */
961	+	#define eat_token(spp,c)        ((next_token(spp)==(c)) ^ !(c) ? *(*(spp))++ : 0)
962	+
963	+	/* Count words from this point in string to '}' */
964	+	static int
965	+	count_values(char *cp)
966	+	{
967	+	int     n = 0;
968	+
969	+	while (next_token(&cp) != '}' && *cp) {
970	+	while (!isspace(*cp) & (*cp != ',') & (*cp != '}'))
971	+	if (!*++cp)
972	+	break;
973	+	++n;
974	+	eat_token(&cp, ',');
975	+	}
976	+	return n;
977	+	}
978	+
979	+	/* Load an array of real numbers, returning total */
980	+	static int
981	+	load_values(char **spp, float *va, int n)
982	+	{
983	+	float   *v = va;
984	+	char    *svnext;
985	+
986	+	while (n-- > 0 && (svnext = fskip(*spp)) != NULL) {
987	+	if ((*v++ = atof(*spp)) < 0)
988	+	v[-1] = 0;
989	+	*spp = svnext;
990	+	eat_token(spp, ',');
991	+	}
992	+	return v - va;
993	+	}
994	+
995	+	/* Load BSDF tree data */
996	+	static SDNode *
997	+	load_tree_data(char **spp, int nd)
998	+	{
999	+	SDNode  *st;
1000	+	int     n;
1001	+
1002	+	if (!eat_token(spp, '{')) {
1003	+	strcpy(SDerrorDetail, "Missing '{' in tensor tree");
1004	+	return NULL;
1005	+	}
1006	+	if (next_token(spp) == '{') {   /* tree branches */
1007	+	st = SDnewNode(nd, -1);
1008	+	if (st == NULL)
1009	+	return NULL;
1010	+	for (n = 0; n < 1<<nd; n++)
1011	+	if ((st->u.t[n] = load_tree_data(spp, nd)) == NULL) {
1012	+	SDfreeTre(st);
1013	+	return NULL;
1014	+	}
1015	+	} else {                        /* else load value grid */
1016	+	int     bsiz;
1017	+	n = count_values(*spp); /* see how big the grid is */
1018	+	for (bsiz = 0; bsiz < 8*sizeof(size_t); bsiz += nd)
1019	+	if (1<<bsiz == n)
1020	+	break;
1021	+	if (bsiz >= 8*sizeof(size_t)) {
1022	+	strcpy(SDerrorDetail, "Illegal value count in tensor tree");
1023	+	return NULL;
1024	+	}
1025	+	st = SDnewNode(nd, bsiz/nd);
1026	+	if (st == NULL)
1027	+	return NULL;
1028	+	if (load_values(spp, st->u.v, n) != n) {
1029	+	strcpy(SDerrorDetail, "Real format error in tensor tree");
1030	+	SDfreeTre(st);
1031	+	return NULL;
1032	+	}
1033	+	}
1034	+	if (!eat_token(spp, '}')) {
1035	+	strcpy(SDerrorDetail, "Missing '}' in tensor tree");
1036	+	SDfreeTre(st);
1037	+	return NULL;
1038	+	}
1039	+	eat_token(spp, ',');
1040	+	return st;
1041	+	}
1042	+
1043	+	/* Compute min. proj. solid angle and max. direct hemispherical scattering */
1044	+	static SDError
1045	+	get_extrema(SDSpectralDF *df)
1046	+	{
1047	+	SDNode  *st = (*(SDTre *)df->comp[0].dist).stc[tt_Y];
1048	+	double  stepWidth, dhemi, bmin[4], bmax[4];
1049	+
1050	+	stepWidth = SDsmallestLeaf(st);
1051	+	if (quantum > stepWidth)        /* adjust quantization factor */
1052	+	quantum = stepWidth;
1053	+	df->minProjSA = M_PI*stepWidth*stepWidth;
1054	+	if (stepWidth < .03125)
1055	+	stepWidth = .03125;     /* 1/32 resolution good enough */
1056	+	df->maxHemi = .0;
1057	+	if (st->ndim == 3) {            /* isotropic BSDF */
1058	+	bmin[1] = bmin[2] = .0;
1059	+	bmax[1] = bmax[2] = 1.;
1060	+	for (bmin[0] = .0; bmin[0] < .5-FTINY; bmin[0] += stepWidth) {
1061	+	bmax[0] = bmin[0] + stepWidth;
1062	+	dhemi = SDavgTreBox(st, bmin, bmax);
1063	+	if (dhemi > df->maxHemi)
1064	+	df->maxHemi = dhemi;
1065	+	}
1066	+	} else if (st->ndim == 4) {     /* anisotropic BSDF */
1067	+	bmin[2] = bmin[3] = .0;
1068	+	bmax[2] = bmax[3] = 1.;
1069	+	for (bmin[0] = .0; bmin[0] < 1.-FTINY; bmin[0] += stepWidth) {
1070	+	bmax[0] = bmin[0] + stepWidth;
1071	+	for (bmin[1] = .0; bmin[1] < 1.-FTINY; bmin[1] += stepWidth) {
1072	+	bmax[1] = bmin[1] + stepWidth;
1073	+	dhemi = SDavgTreBox(st, bmin, bmax);
1074	+	if (dhemi > df->maxHemi)
1075	+	df->maxHemi = dhemi;
1076	+	}
1077	+	}
1078	+	} else
1079	+	return SDEinternal;
1080	+	/* correct hemispherical value */
1081	+	df->maxHemi *= M_PI;
1082	+	return SDEnone;
1083	+	}
1084	+
1085	+	/* Load BSDF distribution for this wavelength */
1086	+	static SDError
1087	+	load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int ndim)
1088	+	{
1089	+	SDSpectralDF    *df;
1090	+	SDTre           *sdt;
1091	+	char            *sdata;
1092	+	/* allocate BSDF component */
1093	+	sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
1094	+	if (!sdata)
1095	+	return SDEnone;
1096	+	/*
1097	+	* Remember that front and back are reversed from WINDOW 6 orientations
1098	+	*/
1099	+	if (!strcasecmp(sdata, "Transmission Front")) {
1100	+	if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(1)) == NULL)
1101	+	return SDEmemory;
1102	+	df = sd->tb;
1103	+	} else if (!strcasecmp(sdata, "Transmission Back")) {
1104	+	if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(1)) == NULL)
1105	+	return SDEmemory;
1106	+	df = sd->tf;
1107	+	} else if (!strcasecmp(sdata, "Reflection Front")) {
1108	+	if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(1)) == NULL)
1109	+	return SDEmemory;
1110	+	df = sd->rb;
1111	+	} else if (!strcasecmp(sdata, "Reflection Back")) {
1112	+	if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(1)) == NULL)
1113	+	return SDEmemory;
1114	+	df = sd->rf;
1115	+	} else
1116	+	return SDEnone;
1117	+	/* get angle bases */
1118	+	sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis"));
1119	+	if (!sdata || strcasecmp(sdata, "LBNL/Shirley-Chiu")) {
1120	+	sprintf(SDerrorDetail, "%s angle basis for BSDF '%s'",
1121	+	!sdata ? "Missing" : "Unsupported", sd->name);
1122	+	return !sdata ? SDEformat : SDEsupport;
1123	+	}
1124	+	if (df->comp[0].dist == NULL) { /* need to allocate BSDF tree? */
1125	+	sdt = (SDTre *)malloc(sizeof(SDTre));
1126	+	if (sdt == NULL)
1127	+	return SDEmemory;
1128	+	if (df == sd->rf)
1129	+	sdt->sidef = SD_FREFL;
1130	+	else if (df == sd->rb)
1131	+	sdt->sidef = SD_BREFL;
1132	+	else if (df == sd->tf)
1133	+	sdt->sidef = SD_FXMIT;
1134	+	else /* df == sd->tb */
1135	+	sdt->sidef = SD_BXMIT;
1136	+	sdt->stc[tt_Y] = sdt->stc[tt_u] = sdt->stc[tt_v] = NULL;
1137	+	df->comp[0].dist = sdt;
1138	+	df->comp[0].func = &SDhandleTre;
1139	+	} else {
1140	+	sdt = (SDTre *)df->comp[0].dist;
1141	+	if (sdt->stc[ct] != NULL) {
1142	+	SDfreeTre(sdt->stc[ct]);
1143	+	sdt->stc[ct] = NULL;
1144	+	}
1145	+	}
1146	+	/* read BSDF data */
1147	+	sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
1148	+	if (!sdata || !next_token(&sdata)) {
1149	+	sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
1150	+	sd->name);
1151	+	return SDEformat;
1152	+	}
1153	+	sdt->stc[ct] = load_tree_data(&sdata, ndim);
1154	+	if (sdt->stc[ct] == NULL)
1155	+	return SDEformat;
1156	+	if (next_token(&sdata)) {       /* check for unconsumed characters */
1157	+	sprintf(SDerrorDetail,
1158	+	"Extra characters at end of ScatteringData in '%s'",
1159	+	sd->name);
1160	+	return SDEformat;
1161	+	}
1162	+	/* flatten branches where possible */
1163	+	sdt->stc[ct] = SDsimplifyTre(sdt->stc[ct]);
1164	+	if (sdt->stc[ct] == NULL)
1165	+	return SDEinternal;
1166	+	/* compute global quantities for Y */
1167	+	return (ct == tt_Y) ? get_extrema(df) : SDEnone;
1168	+	}
1169	+
1170	+	/* Find minimum value in tree */
1171	+	static float
1172	+	SDgetTreMin(const SDNode *st)
1173	+	{
1174	+	float   vmin = FHUGE;
1175	+	int     n;
1176	+
1177	+	if (st->log2GR < 0) {
1178	+	for (n = 1<<st->ndim; n--; ) {
1179	+	float   v = SDgetTreMin(st->u.t[n]);
1180	+	if (v < vmin)
1181	+	vmin = v;
1182	+	}
1183	+	} else {
1184	+	for (n = 1<<(st->ndim*st->log2GR); n--; )
1185	+	if (st->u.v[n] < vmin)
1186	+	vmin = st->u.v[n];
1187	+	}
1188	+	return vmin;
1189	+	}
1190	+
1191	+	/* Subtract the given value from all tree nodes */
1192	+	static void
1193	+	SDsubtractTreVal(SDNode *st, float val)
1194	+	{
1195	+	int     n;
1196	+
1197	+	if (st->log2GR < 0) {
1198	+	for (n = 1<<st->ndim; n--; )
1199	+	SDsubtractTreVal(st->u.t[n], val);
1200	+	} else {
1201	+	for (n = 1<<(st->ndim*st->log2GR); n--; )
1202	+	if ((st->u.v[n] -= val) < 0)
1203	+	st->u.v[n] = .0f;
1204	+	}
1205	+	}
1206	+
1207	+	/* Subtract minimum Y value from BSDF */
1208	+	static double
1209	+	subtract_min_Y(SDNode *st)
1210	+	{
1211	+	const float     vmaxmin = 1.5/M_PI;
1212	+	float           vmin;
1213	+	/* be sure to skip unused portion */
1214	+	if (st->ndim == 3) {
1215	+	int     n;
1216	+	vmin = vmaxmin;
1217	+	if (st->log2GR < 0) {
1218	+	for (n = 0; n < 8; n += 2) {
1219	+	float   v = SDgetTreMin(st->u.t[n]);
1220	+	if (v < vmin)
1221	+	vmin = v;
1222	+	}
1223	+	} else if (st->log2GR) {
1224	+	for (n = 1 << (3*st->log2GR - 1); n--; )
1225	+	if (st->u.v[n] < vmin)
1226	+	vmin = st->u.v[n];
1227	+	} else
1228	+	vmin = st->u.v[0];
1229	+	} else                          /* anisotropic covers entire tree */
1230	+	vmin = SDgetTreMin(st);
1231	+
1232	+	if ((vmin >= vmaxmin) | (vmin <= .01/M_PI))
1233	+	return .0;              /* not worth bothering about */
1234	+
1235	+	SDsubtractTreVal(st, vmin);
1236	+
1237	+	return M_PI * vmin;             /* return hemispherical value */
1238	+	}
1239	+
1240	+	/* Struct used in callback to find RGB extrema */
1241	+	typedef struct {
1242	+	SDNode  **stc;                  /* original Y, u' & v' trees */
1243	+	float   rgb[3];                 /* RGB value */
1244	+	SDNode  *new_stu, *new_stv;     /* replacement u' & v' trees */
1245	+	} SDextRGBs;
1246	+
1247	+	/* Callback to find minimum RGB from Y value plus CIE (u',v') trees */
1248	+	static int
1249	+	get_min_RGB(float yval, const double *cmin, double csiz, void *cptr)
1250	+	{
1251	+	SDextRGBs       *mp = (SDextRGBs *)cptr;
1252	+	double          cmax[SD_MAXDIM];
1253	+	float           rgb[3];
1254	+
1255	+	if (mp->stc[tt_Y]->ndim == 3) {
1256	+	if (cmin[0] + .5*csiz >= .5)
1257	+	return 0;       /* ignore dead half of isotropic */
1258	+	} else
1259	+	cmax[3] = cmin[3] + csiz;
1260	+	cmax[0] = cmin[0] + csiz;
1261	+	cmax[1] = cmin[1] + csiz;
1262	+	cmax[2] = cmin[2] + csiz;
1263	+	/* average RGB color over voxel */
1264	+	SDyuv2rgb(yval, SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1265	+	SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1266	+	/* track smallest components */
1267	+	if (rgb[0] < mp->rgb[0]) mp->rgb[0] = rgb[0];
1268	+	if (rgb[1] < mp->rgb[1]) mp->rgb[1] = rgb[1];
1269	+	if (rgb[2] < mp->rgb[2]) mp->rgb[2] = rgb[2];
1270	+	return 0;
1271	+	}
1272	+
1273	+	/* Callback to build adjusted u' tree */
1274	+	static int
1275	+	adjust_utree(float uprime, const double *cmin, double csiz, void *cptr)
1276	+	{
1277	+	SDextRGBs       *mp = (SDextRGBs *)cptr;
1278	+	double          cmax[SD_MAXDIM];
1279	+	double          yval;
1280	+	float           rgb[3];
1281	+	C_COLOR         clr;
1282	+
1283	+	if (mp->stc[tt_Y]->ndim == 3) {
1284	+	if (cmin[0] + .5*csiz >= .5)
1285	+	return 0;       /* ignore dead half of isotropic */
1286	+	} else
1287	+	cmax[3] = cmin[3] + csiz;
1288	+	cmax[0] = cmin[0] + csiz;
1289	+	cmax[1] = cmin[1] + csiz;
1290	+	cmax[2] = cmin[2] + csiz;
1291	+	/* average RGB color over voxel */
1292	+	SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax), uprime,
1293	+	SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1294	+	/* subtract minimum (& clamp) */
1295	+	if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval;
1296	+	if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval;
1297	+	if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval;
1298	+	c_fromSharpRGB(rgb, &clr);      /* compute new u' for adj. RGB */
1299	+	uprime = 4.*clr.cx/(-2.*clr.cx + 12.*clr.cy + 3.);
1300	+	/* assign in new u' tree */
1301	+	mp->new_stu = SDsetVoxel(mp->new_stu, mp->stc[tt_Y]->ndim,
1302	+	cmin, csiz, uprime);
1303	+	return -(mp->new_stu == NULL);
1304	+	}
1305	+
1306	+	/* Callback to build adjusted v' tree */
1307	+	static int
1308	+	adjust_vtree(float vprime, const double *cmin, double csiz, void *cptr)
1309	+	{
1310	+	SDextRGBs       *mp = (SDextRGBs *)cptr;
1311	+	double          cmax[SD_MAXDIM];
1312	+	double          yval;
1313	+	float           rgb[3];
1314	+	C_COLOR         clr;
1315	+
1316	+	if (mp->stc[tt_Y]->ndim == 3) {
1317	+	if (cmin[0] + .5*csiz >= .5)
1318	+	return 0;       /* ignore dead half of isotropic */
1319	+	} else
1320	+	cmax[3] = cmin[3] + csiz;
1321	+	cmax[0] = cmin[0] + csiz;
1322	+	cmax[1] = cmin[1] + csiz;
1323	+	cmax[2] = cmin[2] + csiz;
1324	+	/* average RGB color over voxel */
1325	+	SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax),
1326	+	SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1327	+	vprime, rgb);
1328	+	/* subtract minimum (& clamp) */
1329	+	if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval;
1330	+	if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval;
1331	+	if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval;
1332	+	c_fromSharpRGB(rgb, &clr);      /* compute new v' for adj. RGB */
1333	+	vprime = 9.*clr.cy/(-2.*clr.cx + 12.*clr.cy + 3.);
1334	+	/* assign in new v' tree */
1335	+	mp->new_stv = SDsetVoxel(mp->new_stv, mp->stc[tt_Y]->ndim,
1336	+	cmin, csiz, vprime);
1337	+	return -(mp->new_stv == NULL);
1338	+	}
1339	+
1340	+	/* Subtract minimum (diffuse) color and return luminance & CIE (x,y) */
1341	+	static double
1342	+	subtract_min_RGB(C_COLOR *cs, SDNode *stc[])
1343	+	{
1344	+	SDextRGBs       my_min;
1345	+	double          ymin;
1346	+
1347	+	my_min.stc = stc;
1348	+	my_min.rgb[0] = my_min.rgb[1] = my_min.rgb[2] = FHUGE;
1349	+	my_min.new_stu = my_min.new_stv = NULL;
1350	+	/* get minimum RGB value */
1351	+	SDtraverseTre(stc[tt_Y], NULL, 0, get_min_RGB, &my_min);
1352	+	/* convert to C_COLOR */
1353	+	ymin =  c_fromSharpRGB(my_min.rgb, cs);
1354	+	if ((ymin >= .5*FHUGE) | (ymin <= .01/M_PI))
1355	+	return .0;              /* close to zero or no tree */
1356	+	/* adjust u' & v' trees */
1357	+	SDtraverseTre(stc[tt_u], NULL, 0, adjust_utree, &my_min);
1358	+	SDtraverseTre(stc[tt_v], NULL, 0, adjust_vtree, &my_min);
1359	+	SDfreeTre(stc[tt_u]); SDfreeTre(stc[tt_v]);
1360	+	stc[tt_u] = SDsimplifyTre(my_min.new_stu);
1361	+	stc[tt_v] = SDsimplifyTre(my_min.new_stv);
1362	+	/* subtract Y & return hemispherical */
1363	+	SDsubtractTreVal(stc[tt_Y], ymin);
1364	+
1365	+	return M_PI * ymin;
1366	+	}
1367	+
1368	+	/* Extract and separate diffuse portion of BSDF */
1369	+	static void
1370	+	extract_diffuse(SDValue *dv, SDSpectralDF *df)
1371	+	{
1372	+	int     n;
1373	+	SDTre   *sdt;
1374	+
1375	+	if (df == NULL || df->ncomp <= 0) {
1376	+	dv->spec = c_dfcolor;
1377	+	dv->cieY = .0;
1378	+	return;
1379	+	}
1380	+	sdt = (SDTre *)df->comp[0].dist;
1381	+	/* subtract minimum color/grayscale */
1382	+	if (sdt->stc[tt_u] != NULL && sdt->stc[tt_v] != NULL) {
1383	+	int     i = 3*(tt_RGB_coef[1] < .001);
1384	+	while (i--) {           /* initialize on first call */
1385	+	float   rgb[3];
1386	+	rgb[0] = rgb[1] = rgb[2] = .0f; rgb[i] = 1.f;
1387	+	tt_RGB_coef[i] = c_fromSharpRGB(rgb, &tt_RGB_prim[i]);
1388	+	}
1389	+	memcpy(df->comp[0].cspec, tt_RGB_prim, sizeof(tt_RGB_prim));
1390	+	dv->cieY = subtract_min_RGB(&dv->spec, sdt->stc);
1391	+	} else {
1392	+	df->comp[0].cspec[0] = dv->spec = c_dfcolor;
1393	+	dv->cieY = subtract_min_Y(sdt->stc[tt_Y]);
1394	+	}
1395	+	df->maxHemi -= dv->cieY;        /* adjust maximum hemispherical */
1396	+
1397	+	c_ccvt(&dv->spec, C_CSXY);      /* make sure (x,y) is set */
1398	+	}
1399	+
1400		/* Load a variable-resolution BSDF tree from an open XML file */
1401		SDError
1402	<	SDloadTre(SDData *sd, ezxml_t fl)
1402	>	SDloadTre(SDData *sd, ezxml_t wtl)
1403		{
1404	<	return SDEsupport;
1404	>	SDError         ec;
1405	>	ezxml_t         wld, wdb;
1406	>	int             rank;
1407	>	char            *txt;
1408	>	/* basic checks and tensor rank */
1409	>	txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
1410	>	"DataDefinition"), "IncidentDataStructure"));
1411	>	if (txt == NULL || !*txt) {
1412	>	sprintf(SDerrorDetail,
1413	>	"BSDF \"%s\": missing IncidentDataStructure",
1414	>	sd->name);
1415	>	return SDEformat;
1416	>	}
1417	>	if (!strcasecmp(txt, "TensorTree3"))
1418	>	rank = 3;
1419	>	else if (!strcasecmp(txt, "TensorTree4"))
1420	>	rank = 4;
1421	>	else {
1422	>	sprintf(SDerrorDetail,
1423	>	"BSDF \"%s\": unsupported IncidentDataStructure",
1424	>	sd->name);
1425	>	return SDEsupport;
1426	>	}
1427	>	/* load BSDF components */
1428	>	for (wld = ezxml_child(wtl, "WavelengthData");
1429	>	wld != NULL; wld = wld->next) {
1430	>	const char      *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
1431	>	int             ct = -1;
1432	>	if (!strcasecmp(cnm, "Visible"))
1433	>	ct = tt_Y;
1434	>	else if (!strcasecmp(cnm, "CIE-u"))
1435	>	ct = tt_u;
1436	>	else if (!strcasecmp(cnm, "CIE-v"))
1437	>	ct = tt_v;
1438	>	else
1439	>	continue;
1440	>	for (wdb = ezxml_child(wld, "WavelengthDataBlock");
1441	>	wdb != NULL; wdb = wdb->next)
1442	>	if ((ec = load_bsdf_data(sd, wdb, ct, rank)) != SDEnone)
1443	>	return ec;
1444	>	}
1445	>	/* separate diffuse components */
1446	>	extract_diffuse(&sd->rLambFront, sd->rf);
1447	>	extract_diffuse(&sd->rLambBack, sd->rb);
1448	>	extract_diffuse(&sd->tLambFront, sd->tf);
1449	>	if (sd->tb != NULL) {
1450	>	extract_diffuse(&sd->tLambBack, sd->tb);
1451	>	if (sd->tf == NULL)
1452	>	sd->tLambFront = sd->tLambBack;
1453	>	} else if (sd->tf != NULL)
1454	>	sd->tLambBack = sd->tLambFront;
1455	>	/* return success */
1456	>	return SDEnone;
1457		}
1458
1459		/* Variable resolution BSDF methods */
1460		const SDFunc SDhandleTre = {
1461	<	NULL,
1462	<	NULL,
1463	<	NULL,
1464	<	NULL,
1465	<	&SDfreeTree,
1461	>	&SDgetTreBSDF,
1462	>	&SDqueryTreProjSA,
1463	>	&SDgetTreCDist,
1464	>	&SDsampTreCDist,
1465	>	&SDFreeBTre,
1466		};

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines