[ViewVC] Diff of: radiance/ray/src/common/bsdf

Comparing ray/src/common/bsdf_t.c (file contents):
Revision 3.8 by greg, Wed Apr 27 20:05:15 2011 UTC vs.
Revision 3.44 by greg, Sun Aug 23 00:56:00 2015 UTC

#	Line 10 \| Line 10 \| static const char RCSid[] = "$Id$";
10		*
11		*/
12
13	+	#define _USE_MATH_DEFINES
14		#include "rtio.h"
15		#include <stdlib.h>
16		#include <math.h>
#	Line 20 \| Line 21 \| static const char RCSid[] = "$Id$";
21		#include "hilbert.h"
22
23		/* Callback function type for SDtraverseTre() */
24	<	typedef int SDtreCallback(float val, const double *cmin,
25	<	double csiz, void *cptr);
25	<
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))-1;
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	<	int nic; /* number of input coordinates */
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 */
#	Line 61 \| Line 73 \| SDnewNode(int nd, int lg)
73		if (lg < 0) {
74		st = (SDNode *)malloc(sizeof(SDNode) +
75		sizeof(st->u.t[0])*((1<<nd) - 1));
76	<	if (st != NULL)
65	<	memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
66	<	} else
67	<	st = (SDNode *)malloc(sizeof(SDNode) +
68	<	sizeof(st->u.v[0])((1 << ndlg) - 1));
69	<
70	<	if (st == NULL) {
71	<	if (lg < 0)
76	>	if (st == NULL) {
77		sprintf(SDerrorDetail,
78		"Cannot allocate %d branch BSDF tree", 1<<nd);
79	<	else
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	>	sizeof(st->u.v[0])((1 << ndlg) - 1));
85	>	if (st == NULL) {
86		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 85 \| Line 97 \| SDnewNode(int nd, int lg)
97		static void
98		SDfreeTre(SDNode *st)
99		{
100	<	int i;
100	>	int n;
101
102		if (st == NULL)
103		return;
104	<	for (i = (st->log2GR < 0) << st->ndim; i--; )
105	<	SDfreeTre(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 */
#	Line 102 \| Line 114 \| *SDFreeBTre(void p)**
114
115		if (sdt == NULL)
116		return;
117	<	SDfreeTre(sdt->st);
117	>	SDfreeTre(sdt->stc[tt_Y]);
118	>	SDfreeTre(sdt->stc[tt_u]);
119	>	SDfreeTre(sdt->stc[tt_v]);
120		free(sdt);
121		}
122
#	Line 125 \| Line 139 \| *fill_grid_branch(float dptr, const float sptr, int n*
139		static float *
140		grid_branch_start(SDNode *st, int n)
141		{
142	<	unsigned skipsiz = 1 << st->log2GR;
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 >> 1;
148	>	vptr += skipsiz;
149		return vptr;
150		}
151
#	Line 151 \| Line 165 \| *SDsimplifyTre(SDNode st)**
165		return NULL; /* propogate error up call stack */
166		match &= (st->u.t[n]->log2GR == st->u.t[0]->log2GR);
167		}
168	<	if (match && st->u.t[0]->log2GR >= 0) {
169	<	SDNode *stn = SDnewNode(st->ndim, st->u.t[0]->log2GR + 1);
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, st->ndim, st->log2GR);
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)
#	Line 187 \| Line 263 \| *SDsmallestLeaf(const SDNode st)**
263		static double
264		SDiterSum(const float va, int nd, int shft, const int imin, const int *imax)
265		{
266	<	const unsigned skipsiz = 1 << nd*shft;
266	>	const unsigned skipsiz = 1 << --nd*shft;
267		double sum = .0;
268		int i;
269	<
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,
200	<	nd-1, shft, 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
#	Line 205 \| Line 282 \| *SDiterSum(const float va, int nd, int shft, const int**
282		static double
283		SDavgTreBox(const SDNode st, const double bmin, const double *bmax)
284		{
208	–	int imin[SD_MAXDIM], imax[SD_MAXDIM];
285		unsigned n;
286		int i;
287
#	Line 215 \| Line 291 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
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 223 \| Line 299 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
299		if (st->log2GR < 0) { /* iterate on subtree */
300		double sum = .0, wsum = 1e-20;
301		double sbmin[SD_MAXDIM], sbmax[SD_MAXDIM], w;
226	–
302		for (n = 1 << st->ndim; n--; ) {
303		w = 1.;
304		for (i = st->ndim; i--; ) {
#	Line 235 \| Line 310 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
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) {
#	Line 243 \| Line 322 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
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 */
248	<	for (i = st->ndim; i--; ) {
249	<	imin[i] = (bmin[i] <= 0) ? 0
250	<	: (int)((1 << st->log2GR)*bmin[i]);
251	<	imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR)
252	<	: (int)((1 << st->log2GR)*bmax[i] + .999999);
253	<	n *= imax[i] - imin[i];
254	<	}
255	<	if (!n)
256	<	return .0;
257	<
258	<	return SDiterSum(st->u.v, st->ndim, st->log2GR, imin, imax) / (double)n;
340	>	return .0;
341		}
342
343		/* Recursive call for SDtraverseTre() */
#	Line 267 \| Line 349 \| *SDdotravTre(const SDNode st, const double pos, int c*
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;
273	–
358		csiz *= .5;
359		for (i = st->ndim; i--; )
360	<	if (1<<i & cmask)
360	>	if (1<<i & cmask) {
361		if (pos[i] < cmin[i] + csiz)
362	<	for (n = 1 << st->ndim; n--; )
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--; )
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;
#	Line 315 \| Line 402 \| *SDdotravTre(const SDNode st, const double pos, int c*
402		clim[i][0] = 0;
403		clim[i][1] = 1 << st->log2GR;
404		}
318	–	/* fill in unused dimensions */
319	–	for (i = SD_MAXDIM; i-- > st->ndim; ) {
320	–	clim[i][0] = 0; clim[i][1] = 1;
321	–	}
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	<	for (cpos[2] = clim[2][0]; cpos[2] < clim[2][1]; cpos[2]++) {
412	<	bmin[3] = cmin[3] + csiz*(cpos[3] = clim[3][0]);
411	>	if (st->ndim == 3) {
412	>	cpos[2] = clim[2][0];
413		n = cpos[0];
414	<	for (i = 1; i < st->ndim; i++)
414	>	for (i = 1; i < 3; i++)
415		n = (n << st->log2GR) + cpos[i];
416	<	for ( ; cpos[3] < clim[3][1]; cpos[3]++) {
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[3] += csiz;
420	>	bmin[2] += csiz;
421		}
422	<	bmin[2] += csiz;
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		}
#	Line 354 \| Line 450 \| static int
450		SDtraverseTre(const SDNode st, const double pos, int cmask,
451		SDtreCallback cf, void cptr)
452		{
357	–	static double czero[SD_MAXDIM];
453		int i;
454		/* check arguments */
455		if ((st == NULL) \| (cf == NULL))
#	Line 379 \| Line 474 \| *SDlookupTre(const SDNode st, const double pos, doubl*
474		hcube[i] = .0;
475		}
476		/* climb the tree */
477	<	while (st->log2GR < 0) {
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--; ) {
#	Line 392 \| Line 487 \| *SDlookupTre(const SDNode st, const double pos, doubl*
487		st = st->u.t[n]; /* avoids tail recursion */
488		pos = spos;
489		}
490	+	if (st == NULL) /* should never happen? */
491	+	return .0;
492		if (st->log2GR == 0) /* short cut */
493		return st->u.v[0];
494		n = t = 0; /* find grid array index */
#	Line 407 \| Line 504 \| *SDlookupTre(const SDNode st, const double pos, doubl*
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.uprimedfact, 4.vprimedfact);
515	+	c_toSharpRGB(&cxy, yval, rgb);
516	+	}
517	+
518		/* Query BSDF value and sample hypercube for the given vectors */
519	<	static float
520	<	SDqueryTre(const SDTre sdt, const FVECT outVec, const FVECT inVec, double hc)
519	>	static int
520	>	SDqueryTre(const SDTre sdt, float coef,
521	>	const FVECT outVec, const FVECT inVec, double *hc)
522		{
523	<	static const FVECT zvec = {.0, .0, 1.};
524	<	FVECT rOutVec;
525	<	double gridPos[4];
523	>	const RREAL *vtmp;
524	>	float yval;
525	>	FVECT rOutVec;
526	>	double gridPos[4];
527
528	+	if (sdt->stc[tt_Y] == NULL) /* paranoia, I hope */
529	+	return 0;
530	+
531		switch (sdt->sidef) { /* whose side are you on? */
532	<	case SD_UFRONT:
532	>	case SD_FREFL:
533		if ((outVec[2] < 0) \| (inVec[2] < 0))
534	<	return -1.;
534	>	return 0;
535		break;
536	<	case SD_UBACK:
536	>	case SD_BREFL:
537		if ((outVec[2] > 0) \| (inVec[2] > 0))
538	<	return -1.;
538	>	return 0;
539		break;
540	<	case SD_XMIT:
541	<	if ((outVec[2] > 0) == (inVec[2] > 0))
542	<	return -1.;
540	>	case SD_FXMIT:
541	>	if (outVec[2] > 0) {
542	>	if (inVec[2] > 0)
543	>	return 0;
544	>	vtmp = outVec; outVec = inVec; inVec = vtmp;
545	>	} else if (inVec[2] < 0)
546	>	return 0;
547		break;
548	+	case SD_BXMIT:
549	+	if (inVec[2] > 0) {
550	+	if (outVec[2] > 0)
551	+	return 0;
552	+	vtmp = outVec; outVec = inVec; inVec = vtmp;
553	+	} else if (outVec[2] < 0)
554	+	return 0;
555	+	break;
556		default:
557	<	return -1.;
557	>	return 0;
558		}
559		/* convert vector coordinates */
560	<	if (sdt->st->ndim == 3) {
561	<	spinvector(rOutVec, outVec, zvec, -atan2(inVec[1],inVec[0]));
562	<	gridPos[0] = .5 - .5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
560	>	if (sdt->stc[tt_Y]->ndim == 3) {
561	>	spinvector(rOutVec, outVec, zvec, -atan2(-inVec[1],-inVec[0]));
562	>	gridPos[0] = (.5-FTINY) -
563	>	.5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
564		SDdisk2square(gridPos+1, rOutVec[0], rOutVec[1]);
565	<	} else if (sdt->st->ndim == 4) {
565	>	} else if (sdt->stc[tt_Y]->ndim == 4) {
566		SDdisk2square(gridPos, -inVec[0], -inVec[1]);
567		SDdisk2square(gridPos+2, outVec[0], outVec[1]);
568		} else
569	<	return -1.; /* should be internal error */
570	<
571	<	return SDlookupTre(sdt->st, gridPos, hc);
569	>	return 0; /* should be internal error */
570	>	/* get BSDF value */
571	>	yval = SDlookupTre(sdt->stc[tt_Y], gridPos, hc);
572	>	if (coef == NULL) /* just getting hypercube? */
573	>	return 1;
574	>	if (sdt->stc[tt_u] == NULL \|\| sdt->stc[tt_v] == NULL) {
575	>	*coef = yval;
576	>	return 1; /* no color */
577	>	}
578	>	/* else decode color */
579	>	SDyuv2rgb(yval, SDlookupTre(sdt->stc[tt_u], gridPos, NULL),
580	>	SDlookupTre(sdt->stc[tt_v], gridPos, NULL), coef);
581	>	coef[0] *= tt_RGB_coef[0];
582	>	coef[1] *= tt_RGB_coef[1];
583	>	coef[2] *= tt_RGB_coef[2];
584	>	return 3;
585		}
586
587		/* Compute non-diffuse component for variable-resolution BSDF */
#	Line 455 \| Line 594 \| SDgetTreBSDF(float coef[SDmaxCh], const FVECT outVec,
594		\|\| sdc->dist == NULL)
595		return 0;
596		/* get nearest BSDF value */
597	<	coef[0] = SDqueryTre((SDTre *)sdc->dist, outVec, inVec, NULL);
459	<	return (coef[0] >= 0); /* monochromatic for now */
597	>	return SDqueryTre((SDTre *)sdc->dist, coef, outVec, inVec, NULL);
598		}
599
600		/* Callback to build cumulative distribution using SDtraverseTre() */
#	Line 465 \| Line 603 \| *build_scaffold(float val, const double cmin, double c**
603		{
604		SDdistScaffold sp = (SDdistScaffold )cptr;
605		int wid = csiz*(double)iwmax + .5;
606	+	double revcmin[2];
607		bitmask_t bmin[2], bmax[2];
608
609	<	cmin += sp->nic; /* skip to output coords */
609	>	if (sp->rev) { /* need to reverse sense? */
610	>	revcmin[0] = 1. - cmin[0] - csiz;
611	>	revcmin[1] = 1. - cmin[1] - csiz;
612	>	cmin = revcmin;
613	>	} else {
614	>	cmin += sp->nic; /* else skip to output coords */
615	>	}
616		if (wid < sp->wmin) /* new minimum width? */
617		sp->wmin = wid;
618		if (wid > sp->wmax) /* new maximum? */
619		sp->wmax = wid;
620		if (sp->alen >= sp->nall) { /* need more space? */
621		struct outdir_s *ndarr;
622	<	sp->nall += 8192;
622	>	sp->nall = (int)(1.5*sp->nall) + 256;
623		ndarr = (struct outdir_s *)realloc(sp->darr,
624		sizeof(struct outdir_s)*sp->nall);
625	<	if (ndarr == NULL)
625	>	if (ndarr == NULL) {
626	>	sprintf(SDerrorDetail,
627	>	"Cannot grow scaffold to %u entries", sp->nall);
628		return -1; /* abort build */
629	+	}
630		sp->darr = ndarr;
631		}
632		/* find Hilbert entry index */
633		bmin[0] = cmin[0]*(double)iwmax + .5;
634		bmin[1] = cmin[1]*(double)iwmax + .5;
635	<	bmax[0] = bmin[0] + wid;
636	<	bmax[1] = bmin[1] + wid;
635	>	bmax[0] = bmin[0] + wid-1;
636	>	bmax[1] = bmin[1] + wid-1;
637		hilbert_box_vtx(2, sizeof(bitmask_t), iwbits, 1, bmin, bmax);
638		sp->darr[sp->alen].hent = hilbert_c2i(2, iwbits, bmin);
639		sp->darr[sp->alen].wid = wid;
#	Line 498 \| Line 646 \| *build_scaffold(float val, const double cmin, double c**
646		static int
647		sscmp(const void p1, const void p2)
648		{
649	<	return (int)(((const struct outdir_s )p1).hent -
650	<	((const struct outdir_s )p2).hent);
649	>	unsigned h1 = ((const struct outdir_s )p1).hent;
650	>	unsigned h2 = ((const struct outdir_s )p2).hent;
651	>
652	>	if (h1 > h2)
653	>	return 1;
654	>	if (h1 < h2)
655	>	return -1;
656	>	return 0;
657		}
658
659		/* Create a new cumulative distribution for the given input direction */
660		static SDTreCDst *
661	<	make_cdist(const SDTre sdt, const double pos)
661	>	make_cdist(const SDTre sdt, const double invec, int rev)
662		{
663		SDdistScaffold myScaffold;
664	+	double pos[4];
665	+	int cmask;
666		SDTreCDst *cd;
667		struct outdir_s *sp;
668		double scale, cursum;
#	Line 514 \| Line 670 \| *make_cdist(const SDTre sdt, const double pos)*
670		/* initialize scaffold */
671		myScaffold.wmin = iwmax;
672		myScaffold.wmax = 0;
673	<	myScaffold.nic = sdt->st->ndim - 2;
673	>	myScaffold.nic = sdt->stc[tt_Y]->ndim - 2;
674	>	myScaffold.rev = rev;
675		myScaffold.alen = 0;
676	<	myScaffold.nall = 8192;
676	>	myScaffold.nall = 512;
677		myScaffold.darr = (struct outdir_s )malloc(sizeof(struct outdir_s)
678		myScaffold.nall);
679		if (myScaffold.darr == NULL)
680		return NULL;
681	+	/* set up traversal */
682	+	cmask = (1<<myScaffold.nic) - 1;
683	+	for (i = myScaffold.nic; i--; )
684	+	pos[i+2*rev] = invec[i];
685	+	cmask <<= 2*rev;
686		/* grow the distribution */
687	<	if (SDtraverseTre(sdt->st, pos, (1<<myScaffold.nic)-1,
688	<	&build_scaffold, &myScaffold) < 0) {
687	>	if (SDtraverseTre(sdt->stc[tt_Y], pos, cmask,
688	>	build_scaffold, &myScaffold) < 0) {
689		free(myScaffold.darr);
690		return NULL;
691		}
#	Line 531 \| Line 693 \| *make_cdist(const SDTre sdt, const double pos)*
693		cd = (SDTreCDst *)malloc(sizeof(SDTreCDst) +
694		sizeof(cd->carr[0])*myScaffold.alen);
695		if (cd == NULL) {
696	+	sprintf(SDerrorDetail,
697	+	"Cannot allocate %u entry cumulative distribution",
698	+	myScaffold.alen);
699		free(myScaffold.darr);
700		return NULL;
701		}
702	+	cd->isodist = (myScaffold.nic == 1);
703		/* sort the distribution */
704		qsort(myScaffold.darr, cd->calen = myScaffold.alen,
705	<	sizeof(struct outdir_s), &sscmp);
705	>	sizeof(struct outdir_s), sscmp);
706
707		/* record input range */
708		scale = myScaffold.wmin / (double)iwmax;
709		for (i = myScaffold.nic; i--; ) {
710	<	cd->clim[i][0] = floor(pos[i]/scale) * scale;
710	>	cd->clim[i][0] = floor(pos[i+2rev]/scale) scale;
711		cd->clim[i][1] = cd->clim[i][0] + scale;
712		}
713	+	if (cd->isodist) { /* avoid issue in SDqueryTreProjSA() */
714	+	cd->clim[1][0] = cd->clim[0][0];
715	+	cd->clim[1][1] = cd->clim[0][1];
716	+	}
717		cd->max_psa = myScaffold.wmax / (double)iwmax;
718		cd->max_psa = cd->max_psa M_PI;
719	<	cd->sidef = sdt->sidef;
719	>	if (rev)
720	>	cd->sidef = (sdt->sidef==SD_BXMIT) ? SD_FXMIT : SD_BXMIT;
721	>	else
722	>	cd->sidef = sdt->sidef;
723		cd->cTotal = 1e-20; /* compute directional total */
724		sp = myScaffold.darr;
725		for (i = myScaffold.alen; i--; sp++)
#	Line 573 \| Line 746 \| *SDgetTreCDist(const FVECT inVec, SDComponent sdc)**
746		{
747		const SDTre *sdt;
748		double inCoord[2];
576	–	int vflags;
749		int i;
750	+	int mode;
751		SDTreCDst cd, cdlast;
752		/* check arguments */
753		if ((inVec == NULL) \| (sdc == NULL) \|\|
754		(sdt = (SDTre *)sdc->dist) == NULL)
755		return NULL;
756	<	if (sdt->st->ndim == 3) /* isotropic BSDF? */
757	<	inCoord[0] = .5 - .5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
758	<	else if (sdt->st->ndim == 4)
759	<	SDdisk2square(inCoord, -inVec[0], -inVec[1]);
760	<	else
756	>	switch (mode = sdt->sidef) { /* check direction */
757	>	case SD_FREFL:
758	>	if (inVec[2] < 0)
759	>	return NULL;
760	>	break;
761	>	case SD_BREFL:
762	>	if (inVec[2] > 0)
763	>	return NULL;
764	>	break;
765	>	case SD_FXMIT:
766	>	if (inVec[2] < 0)
767	>	mode = SD_BXMIT;
768	>	break;
769	>	case SD_BXMIT:
770	>	if (inVec[2] > 0)
771	>	mode = SD_FXMIT;
772	>	break;
773	>	default:
774	>	return NULL;
775	>	}
776	>	if (sdt->stc[tt_Y]->ndim == 3) { /* isotropic BSDF? */
777	>	if (mode != sdt->sidef) /* XXX unhandled reciprocity */
778	>	return &SDemptyCD;
779	>	inCoord[0] = (.5-FTINY) -
780	>	.5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
781	>	} else if (sdt->stc[tt_Y]->ndim == 4) {
782	>	if (mode != sdt->sidef) /* use reciprocity? */
783	>	SDdisk2square(inCoord, inVec[0], inVec[1]);
784	>	else
785	>	SDdisk2square(inCoord, -inVec[0], -inVec[1]);
786	>	} else
787		return NULL; /* should be internal error */
788	+	/* quantize to avoid f.p. errors */
789	+	for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
790	+	inCoord[i] = floor(inCoord[i]/quantum)quantum + .5quantum;
791		cdlast = NULL; /* check for direction in cache list */
792	+	/* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
793		for (cd = (SDTreCDst *)sdc->cdList; cd != NULL;
794	<	cdlast = cd, cd = (SDTreCDst *)cd->next) {
795	<	for (i = sdt->st->ndim - 2; i--; )
794	>	cdlast = cd, cd = cd->next) {
795	>	if (cd->sidef != mode)
796	>	continue;
797	>	for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
798		if ((cd->clim[i][0] > inCoord[i]) \|
799		(inCoord[i] >= cd->clim[i][1]))
800		break;
#	Line 597 \| Line 802 \| *SDgetTreCDist(const FVECT inVec, SDComponent sdc)**
802		break; /* means we have a match */
803		}
804		if (cd == NULL) /* need to create new entry? */
805	<	cdlast = cd = make_cdist(sdt, inCoord);
805	>	cdlast = cd = make_cdist(sdt, inCoord, mode != sdt->sidef);
806		if (cdlast != NULL) { /* move entry to head of cache list */
807		cdlast->next = cd->next;
808	<	cd->next = sdc->cdList;
808	>	cd->next = (SDTreCDst *)sdc->cdList;
809		sdc->cdList = (SDCDst *)cd;
810		}
811	+	/* END MUTEX LOCK */
812		return (SDCDst )cd; / ready to go */
813		}
814
#	Line 619 \| Line 825 \| *SDqueryTreProjSA(double psa, const FVECT v1, const RR**
825		/* get projected solid angle(s) */
826		if (v2 != NULL) {
827		const SDTre sdt = (SDTre )sdc->dist;
828	<	double hcube[SD_MAXDIM];
829	<	if (SDqueryTre(sdt, v1, v2, hcube) < 0) {
828	>	double hcube[SD_MAXDIM+1];
829	>	if (!SDqueryTre(sdt, NULL, v1, v2, hcube)) {
830		strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA");
831		return SDEinternal;
832		}
833	<	myPSA[0] = hcube[sdt->st->ndim];
833	>	myPSA[0] = hcube[sdt->stc[tt_Y]->ndim];
834		myPSA[1] = myPSA[0] = myPSA[0] M_PI;
835		} else {
836		const SDTreCDst cd = (const SDTreCDst )SDgetTreCDist(v1, sdc);
837		if (cd == NULL)
838	<	return SDEmemory;
839	<	myPSA[0] = M_PI * (cd->clim[0][1] - cd->clim[0][0]) *
840	<	(cd->clim[1][1] - cd->clim[1][0]);
841	<	myPSA[1] = cd->max_psa;
838	>	myPSA[0] = myPSA[1] = 0;
839	>	else {
840	>	myPSA[0] = M_PI * (cd->clim[0][1] - cd->clim[0][0]) *
841	>	(cd->clim[1][1] - cd->clim[1][0]);
842	>	myPSA[1] = cd->max_psa;
843	>	}
844		}
845		switch (qflags) { /* record based on flag settings */
846		case SDqueryVal:
#	Line 647 \| Line 855 \| *SDqueryTreProjSA(double psa, const FVECT v1, const RR**
855		psa[1] = myPSA[1];
856		/* fall through */
857		case SDqueryMin:
858	<	if (myPSA[0] < psa[0])
858	>	if ((myPSA[0] > 0) & (myPSA[0] < psa[0]))
859		psa[0] = myPSA[0];
860		break;
861		}
#	Line 663 \| Line 871 \| SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
871		const SDTreCDst cd = (const SDTreCDst )cdp;
872		const unsigned target = randX*cumlmax;
873		bitmask_t hndx, hcoord[2];
874	<	double gpos[3];
874	>	double gpos[3], rotangle;
875		int i, iupper, ilower;
876		/* check arguments */
877		if ((ioVec == NULL) \| (cd == NULL))
878		return SDEargument;
879	+	if (!cd->sidef)
880	+	return SDEnone; /* XXX should never happen */
881		if (ioVec[2] > 0) {
882	<	if (!(cd->sidef & SD_UFRONT))
882	>	if ((cd->sidef != SD_FREFL) & (cd->sidef != SD_FXMIT))
883		return SDEargument;
884	<	} else if (!(cd->sidef & SD_UBACK))
884	>	} else if ((cd->sidef != SD_BREFL) & (cd->sidef != SD_BXMIT))
885		return SDEargument;
886		/* binary search to find position */
887		ilower = 0; iupper = cd->calen;
888		while ((i = (iupper + ilower) >> 1) != ilower)
889	<	if ((long)target >= (long)cd->carr[i].cuml)
889	>	if (target >= cd->carr[i].cuml)
890		ilower = i;
891		else
892		iupper = i;
#	Line 694 \| Line 904 \| SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
904		SDsquare2disk(gpos, gpos[0], gpos[1]);
905		/* compute Z-coordinate */
906		gpos[2] = 1. - gpos[0]gpos[0] - gpos[1]gpos[1];
907	<	if (gpos[2] > 0) /* paranoia, I hope */
698	<	gpos[2] = sqrt(gpos[2]);
907	>	gpos[2] = sqrt(gpos[2]*(gpos[2]>0));
908		/* emit from back? */
909	<	if (ioVec[2] > 0 ^ cd->sidef != SD_XMIT)
909	>	if ((cd->sidef == SD_BREFL) \| (cd->sidef == SD_FXMIT))
910		gpos[2] = -gpos[2];
911	<	VCOPY(ioVec, gpos);
911	>	if (cd->isodist) { /* rotate isotropic sample */
912	>	rotangle = atan2(-ioVec[1],-ioVec[0]);
913	>	spinvector(ioVec, gpos, zvec, rotangle);
914	>	} else
915	>	VCOPY(ioVec, gpos);
916		return SDEnone;
917		}
918
#	Line 712 \| Line 925 \| next_token(char spp)**
925		return **spp;
926		}
927
928	+	/* Advance pointer past matching token (or any token if c==0) */
929	+	#define eat_token(spp,c) ((next_token(spp)==(c)) ^ !(c) ? ((spp))++ : 0)
930	+
931		/* Count words from this point in string to '}' */
932		static int
933		count_values(char *cp)
934		{
935		int n = 0;
936
937	<	while (next_token(&cp) != '}') {
938	<	if (*cp == '{')
939	<	return -1;
940	<	while (cp && !isspace(cp))
725	<	++cp;
937	>	while (next_token(&cp) != '}' && *cp) {
938	>	while (!isspace(cp) & (cp != ',') & (*cp != '}'))
939	>	if (!*++cp)
940	>	break;
941		++n;
942	<	cp += (next_token(&cp) == ',');
942	>	eat_token(&cp, ',');
943		}
944		return n;
945		}
#	Line 737 \| Line 952 \| load_values(char spp, float va, int n)*
952		char *svnext;
953
954		while (n-- > 0 && (svnext = fskip(*spp)) != NULL) {
955	<	v++ = atof(spp);
955	>	if ((v++ = atof(spp)) < 0)
956	>	v[-1] = 0;
957		*spp = svnext;
958	<	*spp += (next_token(spp) == ',');
958	>	eat_token(spp, ',');
959		}
960		return v - va;
961		}
#	Line 751 \| Line 967 \| load_tree_data(char spp, int nd)**
967		SDNode *st;
968		int n;
969
970	<	if (next_token(spp) != '{') {
970	>	if (!eat_token(spp, '{')) {
971		strcpy(SDerrorDetail, "Missing '{' in tensor tree");
972		return NULL;
973		}
758	–	++spp; / in tree, now */
974		if (next_token(spp) == '{') { /* tree branches */
975		st = SDnewNode(nd, -1);
976		if (st == NULL)
#	Line 768 \| Line 983 \| load_tree_data(char spp, int nd)**
983		} else { /* else load value grid */
984		int bsiz;
985		n = count_values(spp); / see how big the grid is */
986	<	if (n <= 0) {
772	<	strcpy(SDerrorDetail, "Bad tensor tree data");
773	<	return NULL;
774	<	}
775	<	for (bsiz = 0; bsiz < 8*sizeof(size_t)-1; bsiz += nd)
986	>	for (bsiz = 0; bsiz < 8*sizeof(size_t); bsiz += nd)
987		if (1<<bsiz == n)
988		break;
989		if (bsiz >= 8*sizeof(size_t)) {
#	Line 788 \| Line 999 \| load_tree_data(char spp, int nd)**
999		return NULL;
1000		}
1001		}
1002	<	if (next_token(spp) != '}') {
1002	>	if (!eat_token(spp, '}')) {
1003		strcpy(SDerrorDetail, "Missing '}' in tensor tree");
1004		SDfreeTre(st);
1005		return NULL;
1006		}
1007	<	++spp; / walk past close and return */
797	<	*spp += (next_token(spp) == ',');
1007	>	eat_token(spp, ',');
1008		return st;
1009		}
1010
#	Line 802 \| Line 1012 \| load_tree_data(char spp, int nd)**
1012		static SDError
1013		get_extrema(SDSpectralDF *df)
1014		{
1015	<	SDNode st = ((SDTre *)df->comp[0].dist).st;
1015	>	SDNode st = ((SDTre *)df->comp[0].dist).stc[tt_Y];
1016		double stepWidth, dhemi, bmin[4], bmax[4];
1017
1018		stepWidth = SDsmallestLeaf(st);
1019	+	if (quantum > stepWidth) /* adjust quantization factor */
1020	+	quantum = stepWidth;
1021		df->minProjSA = M_PIstepWidthstepWidth;
1022		if (stepWidth < .03125)
1023		stepWidth = .03125; /* 1/32 resolution good enough */
#	Line 840 \| Line 1052 \| *get_extrema(SDSpectralDF df)**
1052
1053		/* Load BSDF distribution for this wavelength */
1054		static SDError
1055	<	load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
1055	>	load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int ndim)
1056		{
1057		SDSpectralDF *df;
1058		SDTre *sdt;
1059		char *sdata;
848	–	int i;
1060		/* allocate BSDF component */
1061		sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
1062		if (!sdata)
#	Line 853 \| Line 1064 \| *load_bsdf_data(SDData sd, ezxml_t wdb, int ndim)**
1064		/*
1065		* Remember that front and back are reversed from WINDOW 6 orientations
1066		*/
1067	<	if (!strcasecmp(sdata, "Transmission")) {
1068	<	if (sd->tf != NULL)
858	<	SDfreeSpectralDF(sd->tf);
859	<	if ((sd->tf = SDnewSpectralDF(1)) == NULL)
1067	>	if (!strcasecmp(sdata, "Transmission Front")) {
1068	>	if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(1)) == NULL)
1069		return SDEmemory;
1070	+	df = sd->tb;
1071	+	} else if (!strcasecmp(sdata, "Transmission Back")) {
1072	+	if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(1)) == NULL)
1073	+	return SDEmemory;
1074		df = sd->tf;
1075		} else if (!strcasecmp(sdata, "Reflection Front")) {
1076	<	if (sd->rb != NULL) /* note back-front reversal */
864	<	SDfreeSpectralDF(sd->rb);
865	<	if ((sd->rb = SDnewSpectralDF(1)) == NULL)
1076	>	if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(1)) == NULL)
1077		return SDEmemory;
1078		df = sd->rb;
1079		} else if (!strcasecmp(sdata, "Reflection Back")) {
1080	<	if (sd->rf != NULL) /* note front-back reversal */
870	<	SDfreeSpectralDF(sd->rf);
871	<	if ((sd->rf = SDnewSpectralDF(1)) == NULL)
1080	>	if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(1)) == NULL)
1081		return SDEmemory;
1082		df = sd->rf;
1083		} else
1084		return SDEnone;
876	–	/* XXX should also check "ScatteringDataType" for consistency? */
1085		/* get angle bases */
1086		sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis"));
1087		if (!sdata \|\| strcasecmp(sdata, "LBNL/Shirley-Chiu")) {
#	Line 881 \| Line 1089 \| *load_bsdf_data(SDData sd, ezxml_t wdb, int ndim)**
1089		!sdata ? "Missing" : "Unsupported", sd->name);
1090		return !sdata ? SDEformat : SDEsupport;
1091		}
1092	<	/* allocate BSDF tree */
1093	<	sdt = (SDTre *)malloc(sizeof(SDTre));
1094	<	if (sdt == NULL)
1095	<	return SDEmemory;
1096	<	if (df == sd->rf)
1097	<	sdt->sidef = SD_UFRONT;
1098	<	else if (df == sd->rb)
1099	<	sdt->sidef = SD_UBACK;
1100	<	else
1101	<	sdt->sidef = SD_XMIT;
1102	<	sdt->st = NULL;
1103	<	df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
1104	<	df->comp[0].dist = sdt;
1105	<	df->comp[0].func = &SDhandleTre;
1092	>	if (df->comp[0].dist == NULL) { /* need to allocate BSDF tree? */
1093	>	sdt = (SDTre *)malloc(sizeof(SDTre));
1094	>	if (sdt == NULL)
1095	>	return SDEmemory;
1096	>	if (df == sd->rf)
1097	>	sdt->sidef = SD_FREFL;
1098	>	else if (df == sd->rb)
1099	>	sdt->sidef = SD_BREFL;
1100	>	else if (df == sd->tf)
1101	>	sdt->sidef = SD_FXMIT;
1102	>	else /* df == sd->tb */
1103	>	sdt->sidef = SD_BXMIT;
1104	>	sdt->stc[tt_Y] = sdt->stc[tt_u] = sdt->stc[tt_v] = NULL;
1105	>	df->comp[0].dist = sdt;
1106	>	df->comp[0].func = &SDhandleTre;
1107	>	} else {
1108	>	sdt = (SDTre *)df->comp[0].dist;
1109	>	if (sdt->stc[ct] != NULL) {
1110	>	SDfreeTre(sdt->stc[ct]);
1111	>	sdt->stc[ct] = NULL;
1112	>	}
1113	>	}
1114		/* read BSDF data */
1115		sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
1116		if (!sdata \|\| !next_token(&sdata)) {
#	Line 902 \| Line 1118 \| *load_bsdf_data(SDData sd, ezxml_t wdb, int ndim)**
1118		sd->name);
1119		return SDEformat;
1120		}
1121	<	sdt->st = load_tree_data(&sdata, ndim);
1122	<	if (sdt->st == NULL)
1121	>	sdt->stc[ct] = load_tree_data(&sdata, ndim);
1122	>	if (sdt->stc[ct] == NULL)
1123		return SDEformat;
1124		if (next_token(&sdata)) { /* check for unconsumed characters */
1125		sprintf(SDerrorDetail,
#	Line 912 \| Line 1128 \| *load_bsdf_data(SDData sd, ezxml_t wdb, int ndim)**
1128		return SDEformat;
1129		}
1130		/* flatten branches where possible */
1131	<	sdt->st = SDsimplifyTre(sdt->st);
1132	<	if (sdt->st == NULL)
1131	>	sdt->stc[ct] = SDsimplifyTre(sdt->stc[ct]);
1132	>	if (sdt->stc[ct] == NULL)
1133		return SDEinternal;
1134	<	return get_extrema(df); /* compute global quantities */
1134	>	/* compute global quantities for Y */
1135	>	return (ct == tt_Y) ? get_extrema(df) : SDEnone;
1136		}
1137
1138		/* Find minimum value in tree */
1139		static float
1140		SDgetTreMin(const SDNode *st)
1141		{
1142	<	float vmin = 1./M_PI;
1142	>	float vmin = FHUGE;
1143		int n;
1144
1145		if (st->log2GR < 0) {
#	Line 950 \| Line 1167 \| *SDsubtractTreVal(SDNode st, float val)**
1167		SDsubtractTreVal(st->u.t[n], val);
1168		} else {
1169		for (n = 1<<(st->ndim*st->log2GR); n--; )
1170	<	st->u.v[n] -= val;
1170	>	if ((st->u.v[n] -= val) < 0)
1171	>	st->u.v[n] = .0f;
1172		}
1173		}
1174
1175	<	/* Subtract minimum value from BSDF */
1175	>	/* Subtract minimum Y value from BSDF */
1176		static double
1177	<	subtract_min(SDNode *st)
1177	>	subtract_min_Y(SDNode *st)
1178		{
1179		float vmin;
1180		/* be sure to skip unused portion */
1181	<	if ((st->ndim == 3) & (st->log2GR < 0)) {
1182	<	float v;
965	<	int i;
1181	>	if (st->ndim == 3) {
1182	>	int n;
1183		vmin = 1./M_PI;
1184	<	for (i = 0; i < 4; i++) {
1185	<	v = SDgetTreMin(st->u.t[i]);
1186	<	if (v < vmin)
1187	<	vmin = v;
1188	<	}
1184	>	if (st->log2GR < 0) {
1185	>	for (n = 0; n < 8; n += 2) {
1186	>	float v = SDgetTreMin(st->u.t[n]);
1187	>	if (v < vmin)
1188	>	vmin = v;
1189	>	}
1190	>	} else if (st->log2GR) {
1191	>	for (n = 1 << (3*st->log2GR - 1); n--; )
1192	>	if (st->u.v[n] < vmin)
1193	>	vmin = st->u.v[n];
1194	>	} else
1195	>	vmin = st->u.v[0];
1196		} else /* anisotropic covers entire tree */
1197		vmin = SDgetTreMin(st);
1198
1199	<	if (vmin <= FTINY)
1200	<	return .0;
1199	>	if (vmin <= .01/M_PI)
1200	>	return .0; /* not worth bothering about */
1201
1202		SDsubtractTreVal(st, vmin);
1203
1204		return M_PI * vmin; /* return hemispherical value */
1205		}
1206
1207	+	/* Struct used in callback to find RGB extrema */
1208	+	typedef struct {
1209	+	SDNode *stc; / original Y, u' & v' trees */
1210	+	float rgb[3]; /* RGB value */
1211	+	SDNode new_stu, new_stv; /* replacement u' & v' trees */
1212	+	} SDextRGBs;
1213	+
1214	+	/* Callback to find minimum RGB from Y value plus CIE (u',v') trees */
1215	+	static int
1216	+	get_min_RGB(float yval, const double cmin, double csiz, void cptr)
1217	+	{
1218	+	SDextRGBs mp = (SDextRGBs )cptr;
1219	+	double cmax[SD_MAXDIM];
1220	+	float rgb[3];
1221	+
1222	+	if (mp->stc[tt_Y]->ndim == 3) {
1223	+	if (cmin[0] + .5*csiz >= .5)
1224	+	return 0; /* ignore dead half of isotropic */
1225	+	} else
1226	+	cmax[3] = cmin[3] + csiz;
1227	+	cmax[0] = cmin[0] + csiz;
1228	+	cmax[1] = cmin[1] + csiz;
1229	+	cmax[2] = cmin[2] + csiz;
1230	+	/* average RGB color over voxel */
1231	+	SDyuv2rgb(yval, SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1232	+	SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1233	+	/* track smallest components */
1234	+	if (rgb[0] < mp->rgb[0]) mp->rgb[0] = rgb[0];
1235	+	if (rgb[1] < mp->rgb[1]) mp->rgb[1] = rgb[1];
1236	+	if (rgb[2] < mp->rgb[2]) mp->rgb[2] = rgb[2];
1237	+	return 0;
1238	+	}
1239	+
1240	+	/* Callback to build adjusted u' tree */
1241	+	static int
1242	+	adjust_utree(float uprime, const double cmin, double csiz, void cptr)
1243	+	{
1244	+	SDextRGBs mp = (SDextRGBs )cptr;
1245	+	double cmax[SD_MAXDIM];
1246	+	double yval;
1247	+	float rgb[3];
1248	+	C_COLOR clr;
1249	+
1250	+	if (mp->stc[tt_Y]->ndim == 3) {
1251	+	if (cmin[0] + .5*csiz >= .5)
1252	+	return 0; /* ignore dead half of isotropic */
1253	+	} else
1254	+	cmax[3] = cmin[3] + csiz;
1255	+	cmax[0] = cmin[0] + csiz;
1256	+	cmax[1] = cmin[1] + csiz;
1257	+	cmax[2] = cmin[2] + csiz;
1258	+	/* average RGB color over voxel */
1259	+	SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax), uprime,
1260	+	SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1261	+	/* subtract minimum (& clamp) */
1262	+	if ((rgb[0] -= mp->rgb[0]) < 1e-5yval) rgb[0] = 1e-5yval;
1263	+	if ((rgb[1] -= mp->rgb[1]) < 1e-5yval) rgb[1] = 1e-5yval;
1264	+	if ((rgb[2] -= mp->rgb[2]) < 1e-5yval) rgb[2] = 1e-5yval;
1265	+	c_fromSharpRGB(rgb, &clr); /* compute new u' for adj. RGB */
1266	+	uprime = 4.clr.cx/(-2.clr.cx + 12.*clr.cy + 3.);
1267	+	/* assign in new u' tree */
1268	+	mp->new_stu = SDsetVoxel(mp->new_stu, mp->stc[tt_Y]->ndim,
1269	+	cmin, csiz, uprime);
1270	+	return -(mp->new_stu == NULL);
1271	+	}
1272	+
1273	+	/* Callback to build adjusted v' tree */
1274	+	static int
1275	+	adjust_vtree(float vprime, 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),
1293	+	SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1294	+	vprime, rgb);
1295	+	/* subtract minimum (& clamp) */
1296	+	if ((rgb[0] -= mp->rgb[0]) < 1e-5yval) rgb[0] = 1e-5yval;
1297	+	if ((rgb[1] -= mp->rgb[1]) < 1e-5yval) rgb[1] = 1e-5yval;
1298	+	if ((rgb[2] -= mp->rgb[2]) < 1e-5yval) rgb[2] = 1e-5yval;
1299	+	c_fromSharpRGB(rgb, &clr); /* compute new v' for adj. RGB */
1300	+	vprime = 9.clr.cy/(-2.clr.cx + 12.*clr.cy + 3.);
1301	+	/* assign in new v' tree */
1302	+	mp->new_stv = SDsetVoxel(mp->new_stv, mp->stc[tt_Y]->ndim,
1303	+	cmin, csiz, vprime);
1304	+	return -(mp->new_stv == NULL);
1305	+	}
1306	+
1307	+	/* Subtract minimum (diffuse) color and return luminance & CIE (x,y) */
1308	+	static double
1309	+	subtract_min_RGB(C_COLOR cs, SDNode stc[])
1310	+	{
1311	+	SDextRGBs my_min;
1312	+	double ymin;
1313	+
1314	+	my_min.stc = stc;
1315	+	my_min.rgb[0] = my_min.rgb[1] = my_min.rgb[2] = FHUGE;
1316	+	my_min.new_stu = my_min.new_stv = NULL;
1317	+	/* get minimum RGB value */
1318	+	SDtraverseTre(stc[tt_Y], NULL, 0, get_min_RGB, &my_min);
1319	+	/* convert to C_COLOR */
1320	+	ymin = c_fromSharpRGB(my_min.rgb, cs);
1321	+	if (ymin <= .01/M_PI) /* not worth bothering about? */
1322	+	return .0;
1323	+	/* adjust u' & v' trees */
1324	+	SDtraverseTre(stc[tt_u], NULL, 0, adjust_utree, &my_min);
1325	+	SDtraverseTre(stc[tt_v], NULL, 0, adjust_vtree, &my_min);
1326	+	SDfreeTre(stc[tt_u]); SDfreeTre(stc[tt_v]);
1327	+	stc[tt_u] = SDsimplifyTre(my_min.new_stu);
1328	+	stc[tt_v] = SDsimplifyTre(my_min.new_stv);
1329	+	/* subtract Y & return hemispherical */
1330	+	SDsubtractTreVal(stc[tt_Y], ymin);
1331	+
1332	+	return M_PI * ymin;
1333	+	}
1334	+
1335		/* Extract and separate diffuse portion of BSDF */
1336		static void
1337		extract_diffuse(SDValue dv, SDSpectralDF df)
1338		{
1339		int n;
1340	+	SDTre *sdt;
1341
1342		if (df == NULL \|\| df->ncomp <= 0) {
1343		dv->spec = c_dfcolor;
1344		dv->cieY = .0;
1345		return;
1346		}
1347	<	dv->spec = df->comp[0].cspec[0];
1348	<	dv->cieY = subtract_min(((SDTre )df->comp[n].dist).st);
1349	<	/* in case of multiple components */
1350	<	for (n = df->ncomp; --n; ) {
1351	<	double ymin = subtract_min(((SDTre )df->comp[n].dist).st);
1352	<	c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
1353	<	dv->cieY += ymin;
1347	>	sdt = (SDTre *)df->comp[0].dist;
1348	>	/* subtract minimum color/grayscale */
1349	>	if (sdt->stc[tt_u] != NULL && sdt->stc[tt_v] != NULL) {
1350	>	int i = 3*(tt_RGB_coef[1] < .001);
1351	>	while (i--) { /* initialize on first call */
1352	>	float rgb[3];
1353	>	rgb[0] = rgb[1] = rgb[2] = .0f; rgb[i] = 1.f;
1354	>	tt_RGB_coef[i] = c_fromSharpRGB(rgb, &tt_RGB_prim[i]);
1355	>	}
1356	>	memcpy(df->comp[0].cspec, tt_RGB_prim, sizeof(tt_RGB_prim));
1357	>	dv->cieY = subtract_min_RGB(&dv->spec, sdt->stc);
1358	>	} else {
1359	>	df->comp[0].cspec[0] = dv->spec = c_dfcolor;
1360	>	dv->cieY = subtract_min_Y(sdt->stc[tt_Y]);
1361		}
1362		df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */
1363		/* make sure everything is set */
#	Line 1034 \| Line 1394 \| *SDloadTre(SDData sd, ezxml_t wtl)**
1394		/* load BSDF components */
1395		for (wld = ezxml_child(wtl, "WavelengthData");
1396		wld != NULL; wld = wld->next) {
1397	<	if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
1398	<	"Visible"))
1399	<	continue; /* just visible for now */
1397	>	const char *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
1398	>	int ct = -1;
1399	>	if (!strcasecmp(cnm, "Visible"))
1400	>	ct = tt_Y;
1401	>	else if (!strcasecmp(cnm, "CIE-u"))
1402	>	ct = tt_u;
1403	>	else if (!strcasecmp(cnm, "CIE-v"))
1404	>	ct = tt_v;
1405	>	else
1406	>	continue;
1407		for (wdb = ezxml_child(wld, "WavelengthDataBlock");
1408		wdb != NULL; wdb = wdb->next)
1409	<	if ((ec = load_bsdf_data(sd, wdb, rank)) != SDEnone)
1409	>	if ((ec = load_bsdf_data(sd, wdb, ct, rank)) != SDEnone)
1410		return ec;
1411		}
1412		/* separate diffuse components */
1413		extract_diffuse(&sd->rLambFront, sd->rf);
1414		extract_diffuse(&sd->rLambBack, sd->rb);
1415	<	extract_diffuse(&sd->tLamb, sd->tf);
1415	>	if (sd->tf != NULL)
1416	>	extract_diffuse(&sd->tLamb, sd->tf);
1417	>	if (sd->tb != NULL)
1418	>	extract_diffuse(&sd->tLamb, sd->tb);
1419		/* return success */
1420		return SDEnone;
1421		}

Diff Legend

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

Comparing ray/src/common/bsdf_t.c (file contents): Revision 3.8 by greg, Wed Apr 27 20:05:15 2011 UTC vs. Revision 3.44 by greg, Sun Aug 23 00:56:00 2015 UTC

Diff Legend

Comparing ray/src/common/bsdf_t.c (file contents):
Revision 3.8 by greg, Wed Apr 27 20:05:15 2011 UTC vs.
Revision 3.44 by greg, Sun Aug 23 00:56:00 2015 UTC