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

Comparing ray/src/common/bsdf_t.c (file contents):
Revision 3.6 by greg, Sun Apr 24 19:39:21 2011 UTC vs.
Revision 3.24 by greg, Sun Sep 2 15:33:15 2012 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 24 \| Line 25 \| *typedef int SDtreCallback(float val, const double cmi**
25		double csiz, void *cptr);
26
27		/* reference width maximum (1.0) */
28	<	static const unsigned iwmax = (1<<(sizeof(unsigned)*4))-1;
28	>	static const unsigned iwbits = sizeof(unsigned)*4;
29	>	static const unsigned iwmax = 1<<(sizeof(unsigned)*4);
30	>	/* maximum cumulative value */
31	>	static const unsigned cumlmax = ~0;
32	>	/* constant z-vector */
33	>	static const FVECT zvec = {.0, .0, 1.};
34	>	/* quantization value */
35	>	static double quantum = 1./256.;
36
37		/* Struct used for our distribution-building callback */
38		typedef struct {
39	<	int wmin; /* minimum square size so far */
40	<	int wmax; /* maximum square size */
41	<	int nic; /* number of input coordinates */
42	<	int alen; /* current array length */
43	<	int nall; /* number of allocated entries */
39	>	short nic; /* number of input coordinates */
40	>	short rev; /* reversing query */
41	>	unsigned alen; /* current array length */
42	>	unsigned nall; /* number of allocated entries */
43	>	unsigned wmin; /* minimum square size so far */
44	>	unsigned wmax; /* maximum square size */
45		struct outdir_s {
46		unsigned hent; /* entering Hilbert index */
47		int wid; /* this square size */
#	Line 57 \| Line 66 \| SDnewNode(int nd, int lg)
66		}
67		if (lg < 0) {
68		st = (SDNode *)malloc(sizeof(SDNode) +
69	<	((1<<nd) - 1)*sizeof(st->u.t[0]));
70	<	if (st != NULL)
62	<	memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
63	<	} else
64	<	st = (SDNode *)malloc(sizeof(SDNode) +
65	<	((1 << ndlg) - 1)sizeof(st->u.v[0]));
66	<
67	<	if (st == NULL) {
68	<	if (lg < 0)
69	>	sizeof(st->u.t[0])*((1<<nd) - 1));
70	>	if (st == NULL) {
71		sprintf(SDerrorDetail,
72		"Cannot allocate %d branch BSDF tree", 1<<nd);
73	<	else
73	>	return NULL;
74	>	}
75	>	memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
76	>	} else {
77	>	st = (SDNode *)malloc(sizeof(SDNode) +
78	>	sizeof(st->u.v[0])((1 << ndlg) - 1));
79	>	if (st == NULL) {
80		sprintf(SDerrorDetail,
81		"Cannot allocate %d BSDF leaves", 1 << nd*lg);
82	<	return NULL;
82	>	return NULL;
83	>	}
84		}
85		st->ndim = nd;
86		st->log2GR = lg;
#	Line 82 \| Line 91 \| SDnewNode(int nd, int lg)
91		static void
92		SDfreeTre(SDNode *st)
93		{
94	<	int i;
94	>	int n;
95
96		if (st == NULL)
97		return;
98	<	for (i = (st->log2GR < 0) << st->ndim; i--; )
99	<	SDfreeTre(st->u.t[i]);
100	<	free((void *)st);
98	>	for (n = (st->log2GR < 0) << st->ndim; n--; )
99	>	SDfreeTre(st->u.t[n]);
100	>	free(st);
101		}
102
103		/* Free a variable-resolution BSDF */
#	Line 103 \| Line 112 \| *SDFreeBTre(void p)**
112		free(sdt);
113		}
114
115	+	/* Fill branch's worth of grid values from subtree */
116	+	static void
117	+	fill_grid_branch(float dptr, const float sptr, int nd, int shft)
118	+	{
119	+	unsigned n = 1 << (shft-1);
120	+
121	+	if (!--nd) { /* end on the line */
122	+	memcpy(dptr, sptr, sizeof(dptr)n);
123	+	return;
124	+	}
125	+	while (n--) /* recurse on each slice */
126	+	fill_grid_branch(dptr + (n << shft*nd),
127	+	sptr + (n << (shft-1)*nd), nd, shft);
128	+	}
129	+
130	+	/* Get pointer at appropriate offset for the given branch */
131	+	static float *
132	+	grid_branch_start(SDNode *st, int n)
133	+	{
134	+	unsigned skipsiz = 1 << (st->log2GR - 1);
135	+	float *vptr = st->u.v;
136	+	int i;
137	+
138	+	for (i = st->ndim; i--; skipsiz <<= st->log2GR)
139	+	if (1<<i & n)
140	+	vptr += skipsiz;
141	+	return vptr;
142	+	}
143	+
144	+	/* Simplify (consolidate) a tree by flattening uniform depth regions */
145	+	static SDNode *
146	+	SDsimplifyTre(SDNode *st)
147	+	{
148	+	int match, n;
149	+
150	+	if (st == NULL) /* check for invalid tree */
151	+	return NULL;
152	+	if (st->log2GR >= 0) /* grid just returns unaltered */
153	+	return st;
154	+	match = 1; /* check if grids below match */
155	+	for (n = 0; n < 1<<st->ndim; n++) {
156	+	if ((st->u.t[n] = SDsimplifyTre(st->u.t[n])) == NULL)
157	+	return NULL; /* propogate error up call stack */
158	+	match &= (st->u.t[n]->log2GR == st->u.t[0]->log2GR);
159	+	}
160	+	if (match && (match = st->u.t[0]->log2GR) >= 0) {
161	+	SDNode *stn = SDnewNode(st->ndim, match + 1);
162	+	if (stn == NULL) /* out of memory? */
163	+	return st;
164	+	/* transfer values to new grid */
165	+	for (n = 1 << st->ndim; n--; )
166	+	fill_grid_branch(grid_branch_start(stn, n),
167	+	st->u.t[n]->u.v, stn->ndim, stn->log2GR);
168	+	SDfreeTre(st); /* free old tree */
169	+	st = stn; /* return new one */
170	+	}
171	+	return st;
172	+	}
173	+
174	+	/* Find smallest leaf in tree */
175	+	static double
176	+	SDsmallestLeaf(const SDNode *st)
177	+	{
178	+	if (st->log2GR < 0) { /* tree branches */
179	+	double lmin = 1.;
180	+	int n;
181	+	for (n = 1<<st->ndim; n--; ) {
182	+	double lsiz = SDsmallestLeaf(st->u.t[n]);
183	+	if (lsiz < lmin)
184	+	lmin = lsiz;
185	+	}
186	+	return .5*lmin;
187	+	}
188	+	/* leaf grid width */
189	+	return 1. / (double)(1 << st->log2GR);
190	+	}
191	+
192		/* Add up N-dimensional hypercube array values over the given box */
193		static double
194	<	SDiterSum(const float va, int nd, int siz, const int imin, const int *imax)
194	>	SDiterSum(const float va, int nd, int shft, const int imin, const int *imax)
195		{
196	+	const unsigned skipsiz = 1 << --nd*shft;
197		double sum = .0;
111	–	unsigned skipsiz = 1;
198		int i;
199	<
200	<	for (i = nd; --i > 0; )
201	<	skipsiz *= siz;
199	>
200	>	va += imin skipsiz;
201	>
202		if (skipsiz == 1)
203		for (i = imin; i < imax; i++)
204	<	sum += va[i];
204	>	sum += *va++;
205		else
206	<	for (i = imin; i < imax; i++)
207	<	sum += SDiterSum(va + i*skipsiz,
122	<	nd-1, siz, imin+1, imax+1);
206	>	for (i = imin; i < imax; i++, va += skipsiz)
207	>	sum += SDiterSum(va, nd, shft, imin+1, imax+1);
208		return sum;
209		}
210
#	Line 127 \| Line 212 \| *SDiterSum(const float va, int nd, int siz, const int**
212		static double
213		SDavgTreBox(const SDNode st, const double bmin, const double *bmax)
214		{
130	–	int imin[SD_MAXDIM], imax[SD_MAXDIM];
215		unsigned n;
216		int i;
217
#	Line 137 \| Line 221 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
221		for (i = st->ndim; i--; ) {
222		if (bmin[i] >= 1.)
223		return .0;
224	<	if (bmax[i] <= .0)
224	>	if (bmax[i] <= 0)
225		return .0;
226		if (bmin[i] >= bmax[i])
227		return .0;
#	Line 145 \| Line 229 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
229		if (st->log2GR < 0) { /* iterate on subtree */
230		double sum = .0, wsum = 1e-20;
231		double sbmin[SD_MAXDIM], sbmax[SD_MAXDIM], w;
148	–
232		for (n = 1 << st->ndim; n--; ) {
233		w = 1.;
234		for (i = st->ndim; i--; ) {
#	Line 157 \| Line 240 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
240		}
241		if (sbmin[i] < .0) sbmin[i] = .0;
242		if (sbmax[i] > 1.) sbmax[i] = 1.;
243	+	if (sbmin[i] >= sbmax[i]) {
244	+	w = .0;
245	+	break;
246	+	}
247		w *= sbmax[i] - sbmin[i];
248		}
249		if (w > 1e-10) {
#	Line 165 \| Line 252 \| *SDavgTreBox(const SDNode st, const double bmin, cons*
252		}
253		}
254		return sum / wsum;
255	+	} else { /* iterate over leaves */
256	+	int imin[SD_MAXDIM], imax[SD_MAXDIM];
257	+
258	+	n = 1;
259	+	for (i = st->ndim; i--; ) {
260	+	imin[i] = (bmin[i] <= 0) ? 0 :
261	+	(int)((1 << st->log2GR)*bmin[i]);
262	+	imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR) :
263	+	(int)((1 << st->log2GR)*bmax[i] + .999999);
264	+	n *= imax[i] - imin[i];
265	+	}
266	+	if (n)
267	+	return SDiterSum(st->u.v, st->ndim,
268	+	st->log2GR, imin, imax) / (double)n;
269		}
270	<	n = 1; /* iterate over leaves */
170	<	for (i = st->ndim; i--; ) {
171	<	imin[i] = (bmin[i] <= .0) ? 0
172	<	: (int)((1 << st->log2GR)*bmin[i]);
173	<	imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR)
174	<	: (int)((1 << st->log2GR)*bmax[i] + .999999);
175	<	n *= imax[i] - imin[i];
176	<	}
177	<	if (!n)
178	<	return .0;
179	<
180	<	return SDiterSum(st->u.v, st->ndim, 1 << st->log2GR, imin, imax) /
181	<	(double)n;
270	>	return .0;
271		}
272
273		/* Recursive call for SDtraverseTre() */
#	Line 193 \| Line 282 \| *SDdotravTre(const SDNode st, const double pos, int c*
282		/* in branches? */
283		if (st->log2GR < 0) {
284		unsigned skipmask = 0;
196	–
285		csiz *= .5;
286		for (i = st->ndim; i--; )
287		if (1<<i & cmask)
288		if (pos[i] < cmin[i] + csiz)
289	<	for (n = 1 << st->ndim; n--; )
289	>	for (n = 1 << st->ndim; n--; ) {
290		if (n & 1<<i)
291		skipmask \|= 1<<n;
292	+	}
293		else
294	<	for (n = 1 << st->ndim; n--; )
294	>	for (n = 1 << st->ndim; n--; ) {
295		if (!(n & 1<<i))
296		skipmask \|= 1<<n;
297	+	}
298		for (n = 1 << st->ndim; n--; ) {
299		if (1<<n & skipmask)
300		continue;
#	Line 238 \| Line 328 \| *SDdotravTre(const SDNode st, const double pos, int c*
328		clim[i][0] = 0;
329		clim[i][1] = 1 << st->log2GR;
330		}
241	–	/* fill in unused dimensions */
242	–	for (i = SD_MAXDIM; i-- > st->ndim; ) {
243	–	clim[i][0] = 0; clim[i][1] = 1;
244	–	}
331		#if (SD_MAXDIM == 4)
332		bmin[0] = cmin[0] + csiz*clim[0][0];
333		for (cpos[0] = clim[0][0]; cpos[0] < clim[0][1]; cpos[0]++) {
334		bmin[1] = cmin[1] + csiz*clim[1][0];
335		for (cpos[1] = clim[1][0]; cpos[1] < clim[1][1]; cpos[1]++) {
336		bmin[2] = cmin[2] + csiz*clim[2][0];
337	<	for (cpos[2] = clim[2][0]; cpos[2] < clim[2][1]; cpos[2]++) {
338	<	bmin[3] = cmin[3] + csiz*(cpos[3] = clim[3][0]);
337	>	if (st->ndim == 3) {
338	>	cpos[2] = clim[2][0];
339		n = cpos[0];
340	<	for (i = 1; i < st->ndim; i++)
340	>	for (i = 1; i < 3; i++)
341		n = (n << st->log2GR) + cpos[i];
342	<	for ( ; cpos[3] < clim[3][1]; cpos[3]++) {
342	>	for ( ; cpos[2] < clim[2][1]; cpos[2]++) {
343		rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr);
344		if (rv < 0)
345		return rv;
346	<	bmin[3] += csiz;
346	>	bmin[2] += csiz;
347		}
348	<	bmin[2] += csiz;
348	>	} else {
349	>	for (cpos[2] = clim[2][0]; cpos[2] < clim[2][1]; cpos[2]++) {
350	>	bmin[3] = cmin[3] + csiz*(cpos[3] = clim[3][0]);
351	>	n = cpos[0];
352	>	for (i = 1; i < 4; i++)
353	>	n = (n << st->log2GR) + cpos[i];
354	>	for ( ; cpos[3] < clim[3][1]; cpos[3]++) {
355	>	rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr);
356	>	if (rv < 0)
357	>	return rv;
358	>	bmin[3] += csiz;
359	>	}
360	>	bmin[2] += csiz;
361	>	}
362		}
363		bmin[1] += csiz;
364		}
#	Line 334 \| Line 433 \| *SDlookupTre(const SDNode st, const double pos, doubl*
433		static float
434		SDqueryTre(const SDTre sdt, const FVECT outVec, const FVECT inVec, double hc)
435		{
436	<	static const FVECT zvec = {.0, .0, 1.};
437	<	FVECT rOutVec;
438	<	double gridPos[4];
439	<	/* check transmission */
440	<	if (!sdt->isxmit ^ outVec[2] > 0 ^ inVec[2] > 0)
436	>	const RREAL *vtmp;
437	>	FVECT rOutVec;
438	>	double gridPos[4];
439	>
440	>	switch (sdt->sidef) { /* whose side are you on? */
441	>	case SD_FREFL:
442	>	if ((outVec[2] < 0) \| (inVec[2] < 0))
443	>	return -1.;
444	>	break;
445	>	case SD_BREFL:
446	>	if ((outVec[2] > 0) \| (inVec[2] > 0))
447	>	return -1.;
448	>	break;
449	>	case SD_FXMIT:
450	>	if (outVec[2] > 0) {
451	>	if (inVec[2] > 0)
452	>	return -1.;
453	>	vtmp = outVec; outVec = inVec; inVec = vtmp;
454	>	} else if (inVec[2] < 0)
455	>	return -1.;
456	>	break;
457	>	case SD_BXMIT:
458	>	if (inVec[2] > 0) {
459	>	if (outVec[2] > 0)
460	>	return -1.;
461	>	vtmp = outVec; outVec = inVec; inVec = vtmp;
462	>	} else if (outVec[2] < 0)
463	>	return -1.;
464	>	break;
465	>	default:
466		return -1.;
467	+	}
468		/* convert vector coordinates */
469		if (sdt->st->ndim == 3) {
470	<	spinvector(rOutVec, outVec, zvec, -atan2(inVec[1],inVec[0]));
470	>	spinvector(rOutVec, outVec, zvec, -atan2(-inVec[1],-inVec[0]));
471		gridPos[0] = .5 - .5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
472		SDdisk2square(gridPos+1, rOutVec[0], rOutVec[1]);
473		} else if (sdt->st->ndim == 4) {
#	Line 376 \| Line 501 \| *build_scaffold(float val, const double cmin, double c**
501		int wid = csiz*(double)iwmax + .5;
502		bitmask_t bmin[2], bmax[2];
503
504	<	cmin += sp->nic; /* skip to output coords */
504	>	cmin += sp->nic * !sp->rev; /* skip to output coords */
505		if (wid < sp->wmin) /* new minimum width? */
506		sp->wmin = wid;
507		if (wid > sp->wmax) /* new maximum? */
508		sp->wmax = wid;
509		if (sp->alen >= sp->nall) { /* need more space? */
510		struct outdir_s *ndarr;
511	<	sp->nall += 8192;
511	>	sp->nall += 1024;
512		ndarr = (struct outdir_s *)realloc(sp->darr,
513		sizeof(struct outdir_s)*sp->nall);
514	<	if (ndarr == NULL)
514	>	if (ndarr == NULL) {
515	>	sprintf(SDerrorDetail,
516	>	"Cannot grow scaffold to %u entries", sp->nall);
517		return -1; /* abort build */
518	+	}
519		sp->darr = ndarr;
520		}
521		/* find Hilbert entry index */
522		bmin[0] = cmin[0]*(double)iwmax + .5;
523		bmin[1] = cmin[1]*(double)iwmax + .5;
524	<	bmax[0] = bmin[0] + wid;
525	<	bmax[1] = bmin[1] + wid;
526	<	hilbert_box_vtx(2, sizeof(bitmask_t), sizeof(unsigned)*4,
527	<	1, bmin, bmax);
400	<	sp->darr[sp->alen].hent = hilbert_c2i(2, sizeof(unsigned)*4, bmin);
524	>	bmax[0] = bmin[0] + wid-1;
525	>	bmax[1] = bmin[1] + wid-1;
526	>	hilbert_box_vtx(2, sizeof(bitmask_t), iwbits, 1, bmin, bmax);
527	>	sp->darr[sp->alen].hent = hilbert_c2i(2, iwbits, bmin);
528		sp->darr[sp->alen].wid = wid;
529		sp->darr[sp->alen].bsdf = val;
530		sp->alen++; /* on to the next entry */
#	Line 408 \| Line 535 \| *build_scaffold(float val, const double cmin, double c**
535		static int
536		sscmp(const void p1, const void p2)
537		{
538	<	return (int)(((const struct outdir_s )p1).hent -
539	<	((const struct outdir_s )p2).hent);
538	>	unsigned h1 = ((const struct outdir_s )p1).hent;
539	>	unsigned h2 = ((const struct outdir_s )p2).hent;
540	>
541	>	if (h1 > h2)
542	>	return 1;
543	>	if (h1 < h2)
544	>	return -1;
545	>	return 0;
546		}
547
548		/* Create a new cumulative distribution for the given input direction */
549		static SDTreCDst *
550	<	make_cdist(const SDTre sdt, const double pos)
550	>	make_cdist(const SDTre sdt, const double invec, int rev)
551		{
419	–	const unsigned cumlmax = ~0;
552		SDdistScaffold myScaffold;
553	+	double pos[4];
554	+	int cmask;
555		SDTreCDst *cd;
556		struct outdir_s *sp;
557		double scale, cursum;
#	Line 426 \| Line 560 \| *make_cdist(const SDTre sdt, const double pos)*
560		myScaffold.wmin = iwmax;
561		myScaffold.wmax = 0;
562		myScaffold.nic = sdt->st->ndim - 2;
563	+	myScaffold.rev = rev;
564		myScaffold.alen = 0;
565	<	myScaffold.nall = 8192;
565	>	myScaffold.nall = 512;
566		myScaffold.darr = (struct outdir_s )malloc(sizeof(struct outdir_s)
567		myScaffold.nall);
568		if (myScaffold.darr == NULL)
569		return NULL;
570	+	/* set up traversal */
571	+	cmask = (1<<myScaffold.nic) - 1;
572	+	for (i = myScaffold.nic; i--; )
573	+	pos[i+2*rev] = invec[i];
574	+	cmask <<= 2*rev;
575		/* grow the distribution */
576	<	if (SDtraverseTre(sdt->st, pos, (1<<myScaffold.nic)-1,
576	>	if (SDtraverseTre(sdt->st, pos, cmask,
577		&build_scaffold, &myScaffold) < 0) {
578		free(myScaffold.darr);
579		return NULL;
#	Line 442 \| Line 582 \| *make_cdist(const SDTre sdt, const double pos)*
582		cd = (SDTreCDst *)malloc(sizeof(SDTreCDst) +
583		sizeof(cd->carr[0])*myScaffold.alen);
584		if (cd == NULL) {
585	+	sprintf(SDerrorDetail,
586	+	"Cannot allocate %u entry cumulative distribution",
587	+	myScaffold.alen);
588		free(myScaffold.darr);
589		return NULL;
590		}
591	+	cd->isodist = (myScaffold.nic == 1);
592		/* sort the distribution */
593		qsort(myScaffold.darr, cd->calen = myScaffold.alen,
594		sizeof(struct outdir_s), &sscmp);
595
596		/* record input range */
597	<	scale = (double)myScaffold.wmin / iwmax;
597	>	scale = myScaffold.wmin / (double)iwmax;
598		for (i = myScaffold.nic; i--; ) {
599	<	cd->clim[i][0] = floor(pos[i]/scale + .5) * scale;
599	>	cd->clim[i][0] = floor(pos[i]/scale) * scale;
600		cd->clim[i][1] = cd->clim[i][0] + scale;
601		}
602	+	if (cd->isodist) { /* avoid issue in SDqueryTreProjSA() */
603	+	cd->clim[1][0] = cd->clim[0][0];
604	+	cd->clim[1][1] = cd->clim[0][1];
605	+	}
606		cd->max_psa = myScaffold.wmax / (double)iwmax;
607		cd->max_psa = cd->max_psa M_PI;
608	<	cd->isxmit = sdt->isxmit;
608	>	if (rev)
609	>	cd->sidef = (sdt->sidef==SD_BXMIT) ? SD_FXMIT : SD_BXMIT;
610	>	else
611	>	cd->sidef = sdt->sidef;
612		cd->cTotal = 1e-20; /* compute directional total */
613		sp = myScaffold.darr;
614		for (i = myScaffold.alen; i--; sp++)
#	Line 466 \| Line 617 \| *make_cdist(const SDTre sdt, const double pos)*
617		scale = (double)cumlmax / cd->cTotal;
618		sp = myScaffold.darr;
619		for (i = 0; i < cd->calen; i++, sp++) {
620	+	cd->carr[i].hndx = sp->hent;
621		cd->carr[i].cuml = scale*cursum + .5;
622		cursum += sp->bsdf * (double)sp->wid * sp->wid;
623		}
#	Line 483 \| Line 635 \| *SDgetTreCDist(const FVECT inVec, SDComponent sdc)**
635		{
636		const SDTre *sdt;
637		double inCoord[2];
486	–	int vflags;
638		int i;
639	+	int mode;
640		SDTreCDst cd, cdlast;
641		/* check arguments */
642		if ((inVec == NULL) \| (sdc == NULL) \|\|
643		(sdt = (SDTre *)sdc->dist) == NULL)
644		return NULL;
645	<	if (sdt->st->ndim == 3) /* isotropic BSDF? */
645	>	switch (mode = sdt->sidef) { /* check direction */
646	>	case SD_FREFL:
647	>	if (inVec[2] < 0)
648	>	return NULL;
649	>	break;
650	>	case SD_BREFL:
651	>	if (inVec[2] > 0)
652	>	return NULL;
653	>	break;
654	>	case SD_FXMIT:
655	>	if (inVec[2] < 0)
656	>	mode = SD_BXMIT;
657	>	break;
658	>	case SD_BXMIT:
659	>	if (inVec[2] > 0)
660	>	mode = SD_FXMIT;
661	>	break;
662	>	default:
663	>	return NULL;
664	>	}
665	>	if (sdt->st->ndim == 3) { /* isotropic BSDF? */
666	>	if (mode != sdt->sidef) /* XXX unhandled reciprocity */
667	>	return &SDemptyCD;
668		inCoord[0] = .5 - .5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
669	<	else if (sdt->st->ndim == 4)
669	>	} else if (sdt->st->ndim == 4) {
670		SDdisk2square(inCoord, -inVec[0], -inVec[1]);
671	<	else
671	>	} else
672		return NULL; /* should be internal error */
673	+	/* quantize to avoid f.p. errors */
674	+	for (i = sdt->st->ndim - 2; i--; )
675	+	inCoord[i] = floor(inCoord[i]/quantum)quantum + .5quantum;
676		cdlast = NULL; /* check for direction in cache list */
677		for (cd = (SDTreCDst *)sdc->cdList; cd != NULL;
678	<	cdlast = cd, cd = (SDTreCDst *)cd->next) {
678	>	cdlast = cd, cd = cd->next) {
679	>	if (cd->sidef != mode)
680	>	continue;
681		for (i = sdt->st->ndim - 2; i--; )
682		if ((cd->clim[i][0] > inCoord[i]) \|
683		(inCoord[i] >= cd->clim[i][1]))
#	Line 507 \| Line 686 \| *SDgetTreCDist(const FVECT inVec, SDComponent sdc)**
686		break; /* means we have a match */
687		}
688		if (cd == NULL) /* need to create new entry? */
689	<	cdlast = cd = make_cdist(sdt, inCoord);
689	>	cdlast = cd = make_cdist(sdt, inCoord, mode != sdt->sidef);
690		if (cdlast != NULL) { /* move entry to head of cache list */
691		cdlast->next = cd->next;
692	<	cd->next = sdc->cdList;
692	>	cd->next = (SDTreCDst *)sdc->cdList;
693		sdc->cdList = (SDCDst *)cd;
694		}
695		return (SDCDst )cd; / ready to go */
#	Line 531 \| Line 710 \| *SDqueryTreProjSA(double psa, const FVECT v1, const RR**
710		const SDTre sdt = (SDTre )sdc->dist;
711		double hcube[SD_MAXDIM];
712		if (SDqueryTre(sdt, v1, v2, hcube) < 0) {
713	<	if (qflags == SDqueryVal)
714	<	*psa = M_PI;
536	<	return SDEnone;
713	>	strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA");
714	>	return SDEinternal;
715		}
716		myPSA[0] = hcube[sdt->st->ndim];
717		myPSA[1] = myPSA[0] = myPSA[0] M_PI;
#	Line 571 \| Line 749 \| SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
749		{
750		const unsigned nBitsC = 4*sizeof(bitmask_t);
751		const unsigned nExtraBits = 8*(sizeof(bitmask_t)-sizeof(unsigned));
574	–	const unsigned maxval = ~0;
752		const SDTreCDst cd = (const SDTreCDst )cdp;
753	<	const unsigned target = randX*maxval;
753	>	const unsigned target = randX*cumlmax;
754		bitmask_t hndx, hcoord[2];
755	<	double gpos[3];
755	>	double gpos[3], rotangle;
756		int i, iupper, ilower;
757		/* check arguments */
758		if ((ioVec == NULL) \| (cd == NULL))
759		return SDEargument;
760	+	if (!cd->sidef)
761	+	return SDEnone; /* XXX should never happen */
762	+	if (ioVec[2] > 0) {
763	+	if ((cd->sidef != SD_FREFL) & (cd->sidef != SD_FXMIT))
764	+	return SDEargument;
765	+	} else if ((cd->sidef != SD_BREFL) & (cd->sidef != SD_BXMIT))
766	+	return SDEargument;
767		/* binary search to find position */
768		ilower = 0; iupper = cd->calen;
769		while ((i = (iupper + ilower) >> 1) != ilower)
770	<	if ((long)target >= (long)cd->carr[i].cuml)
770	>	if (target >= cd->carr[i].cuml)
771		ilower = i;
772		else
773		iupper = i;
774		/* localize random position */
775	<	randX = (randX*maxval - cd->carr[ilower].cuml) /
775	>	randX = (randX*cumlmax - cd->carr[ilower].cuml) /
776		(double)(cd->carr[iupper].cuml - cd->carr[ilower].cuml);
777		/* index in longer Hilbert curve */
778		hndx = (randXcd->carr[iupper].hndx + (1.-randX)cd->carr[ilower].hndx)
#	Line 599 \| Line 783 \| SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
783		gpos[i] = ((double)hcoord[i] + rand()*(1./(RAND_MAX+.5))) /
784		(double)((bitmask_t)1 << nBitsC);
785		SDsquare2disk(gpos, gpos[0], gpos[1]);
786	+	/* compute Z-coordinate */
787		gpos[2] = 1. - gpos[0]gpos[0] - gpos[1]gpos[1];
788		if (gpos[2] > 0) /* paranoia, I hope */
789		gpos[2] = sqrt(gpos[2]);
790	<	if (ioVec[2] > 0 ^ !cd->isxmit)
790	>	/* emit from back? */
791	>	if ((cd->sidef == SD_BREFL) \| (cd->sidef == SD_FXMIT))
792		gpos[2] = -gpos[2];
793	<	VCOPY(ioVec, gpos);
793	>	if (cd->isodist) { /* rotate isotropic result */
794	>	rotangle = atan2(-ioVec[1],-ioVec[0]);
795	>	VCOPY(ioVec, gpos);
796	>	spinvector(ioVec, ioVec, zvec, rotangle);
797	>	} else
798	>	VCOPY(ioVec, gpos);
799		return SDEnone;
800		}
801
802	+	/* Advance pointer to the next non-white character in the string (or nul) */
803	+	static int
804	+	next_token(char **spp)
805	+	{
806	+	while (isspace(**spp))
807	+	++*spp;
808	+	return **spp;
809	+	}
810	+
811	+	/* Advance pointer past matching token (or any token if c==0) */
812	+	#define eat_token(spp,c) (next_token(spp)==(c) ^ !(c) ? ((spp))++ : 0)
813	+
814	+	/* Count words from this point in string to '}' */
815	+	static int
816	+	count_values(char *cp)
817	+	{
818	+	int n = 0;
819	+
820	+	while (next_token(&cp) != '}' && *cp) {
821	+	while (!isspace(cp) & (cp != ',') & (*cp != '}'))
822	+	if (!*++cp)
823	+	break;
824	+	++n;
825	+	eat_token(&cp, ',');
826	+	}
827	+	return n;
828	+	}
829	+
830	+	/* Load an array of real numbers, returning total */
831	+	static int
832	+	load_values(char *spp, float va, int n)
833	+	{
834	+	float *v = va;
835	+	char *svnext;
836	+
837	+	while (n-- > 0 && (svnext = fskip(*spp)) != NULL) {
838	+	v++ = atof(spp);
839	+	*spp = svnext;
840	+	eat_token(spp, ',');
841	+	}
842	+	return v - va;
843	+	}
844	+
845	+	/* Load BSDF tree data */
846	+	static SDNode *
847	+	load_tree_data(char **spp, int nd)
848	+	{
849	+	SDNode *st;
850	+	int n;
851	+
852	+	if (!eat_token(spp, '{')) {
853	+	strcpy(SDerrorDetail, "Missing '{' in tensor tree");
854	+	return NULL;
855	+	}
856	+	if (next_token(spp) == '{') { /* tree branches */
857	+	st = SDnewNode(nd, -1);
858	+	if (st == NULL)
859	+	return NULL;
860	+	for (n = 0; n < 1<<nd; n++)
861	+	if ((st->u.t[n] = load_tree_data(spp, nd)) == NULL) {
862	+	SDfreeTre(st);
863	+	return NULL;
864	+	}
865	+	} else { /* else load value grid */
866	+	int bsiz;
867	+	n = count_values(spp); / see how big the grid is */
868	+	for (bsiz = 0; bsiz < 8*sizeof(size_t); bsiz += nd)
869	+	if (1<<bsiz == n)
870	+	break;
871	+	if (bsiz >= 8*sizeof(size_t)) {
872	+	strcpy(SDerrorDetail, "Illegal value count in tensor tree");
873	+	return NULL;
874	+	}
875	+	st = SDnewNode(nd, bsiz/nd);
876	+	if (st == NULL)
877	+	return NULL;
878	+	if (load_values(spp, st->u.v, n) != n) {
879	+	strcpy(SDerrorDetail, "Real format error in tensor tree");
880	+	SDfreeTre(st);
881	+	return NULL;
882	+	}
883	+	}
884	+	if (!eat_token(spp, '}')) {
885	+	strcpy(SDerrorDetail, "Missing '}' in tensor tree");
886	+	SDfreeTre(st);
887	+	return NULL;
888	+	}
889	+	eat_token(spp, ',');
890	+	return st;
891	+	}
892	+
893	+	/* Compute min. proj. solid angle and max. direct hemispherical scattering */
894	+	static SDError
895	+	get_extrema(SDSpectralDF *df)
896	+	{
897	+	SDNode st = ((SDTre *)df->comp[0].dist).st;
898	+	double stepWidth, dhemi, bmin[4], bmax[4];
899	+
900	+	stepWidth = SDsmallestLeaf(st);
901	+	if (quantum > stepWidth) /* adjust quantization factor */
902	+	quantum = stepWidth;
903	+	df->minProjSA = M_PIstepWidthstepWidth;
904	+	if (stepWidth < .03125)
905	+	stepWidth = .03125; /* 1/32 resolution good enough */
906	+	df->maxHemi = .0;
907	+	if (st->ndim == 3) { /* isotropic BSDF */
908	+	bmin[1] = bmin[2] = .0;
909	+	bmax[1] = bmax[2] = 1.;
910	+	for (bmin[0] = .0; bmin[0] < .5-FTINY; bmin[0] += stepWidth) {
911	+	bmax[0] = bmin[0] + stepWidth;
912	+	dhemi = SDavgTreBox(st, bmin, bmax);
913	+	if (dhemi > df->maxHemi)
914	+	df->maxHemi = dhemi;
915	+	}
916	+	} else if (st->ndim == 4) { /* anisotropic BSDF */
917	+	bmin[2] = bmin[3] = .0;
918	+	bmax[2] = bmax[3] = 1.;
919	+	for (bmin[0] = .0; bmin[0] < 1.-FTINY; bmin[0] += stepWidth) {
920	+	bmax[0] = bmin[0] + stepWidth;
921	+	for (bmin[1] = .0; bmin[1] < 1.-FTINY; bmin[1] += stepWidth) {
922	+	bmax[1] = bmin[1] + stepWidth;
923	+	dhemi = SDavgTreBox(st, bmin, bmax);
924	+	if (dhemi > df->maxHemi)
925	+	df->maxHemi = dhemi;
926	+	}
927	+	}
928	+	} else
929	+	return SDEinternal;
930	+	/* correct hemispherical value */
931	+	df->maxHemi *= M_PI;
932	+	return SDEnone;
933	+	}
934	+
935	+	/* Load BSDF distribution for this wavelength */
936	+	static SDError
937	+	load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
938	+	{
939	+	SDSpectralDF *df;
940	+	SDTre *sdt;
941	+	char *sdata;
942	+	/* allocate BSDF component */
943	+	sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
944	+	if (!sdata)
945	+	return SDEnone;
946	+	/*
947	+	* Remember that front and back are reversed from WINDOW 6 orientations
948	+	*/
949	+	if (!strcasecmp(sdata, "Transmission Front")) {
950	+	if (sd->tf != NULL)
951	+	SDfreeSpectralDF(sd->tf);
952	+	if ((sd->tf = SDnewSpectralDF(1)) == NULL)
953	+	return SDEmemory;
954	+	df = sd->tf;
955	+	} else if (!strcasecmp(sdata, "Transmission Back")) {
956	+	if (sd->tb != NULL)
957	+	SDfreeSpectralDF(sd->tb);
958	+	if ((sd->tb = SDnewSpectralDF(1)) == NULL)
959	+	return SDEmemory;
960	+	df = sd->tb;
961	+	} else if (!strcasecmp(sdata, "Reflection Front")) {
962	+	if (sd->rb != NULL) /* note back-front reversal */
963	+	SDfreeSpectralDF(sd->rb);
964	+	if ((sd->rb = SDnewSpectralDF(1)) == NULL)
965	+	return SDEmemory;
966	+	df = sd->rb;
967	+	} else if (!strcasecmp(sdata, "Reflection Back")) {
968	+	if (sd->rf != NULL) /* note front-back reversal */
969	+	SDfreeSpectralDF(sd->rf);
970	+	if ((sd->rf = SDnewSpectralDF(1)) == NULL)
971	+	return SDEmemory;
972	+	df = sd->rf;
973	+	} else
974	+	return SDEnone;
975	+	/* XXX should also check "ScatteringDataType" for consistency? */
976	+	/* get angle bases */
977	+	sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis"));
978	+	if (!sdata \|\| strcasecmp(sdata, "LBNL/Shirley-Chiu")) {
979	+	sprintf(SDerrorDetail, "%s angle basis for BSDF '%s'",
980	+	!sdata ? "Missing" : "Unsupported", sd->name);
981	+	return !sdata ? SDEformat : SDEsupport;
982	+	}
983	+	/* allocate BSDF tree */
984	+	sdt = (SDTre *)malloc(sizeof(SDTre));
985	+	if (sdt == NULL)
986	+	return SDEmemory;
987	+	if (df == sd->rf)
988	+	sdt->sidef = SD_FREFL;
989	+	else if (df == sd->rb)
990	+	sdt->sidef = SD_BREFL;
991	+	else if (df == sd->tf)
992	+	sdt->sidef = SD_FXMIT;
993	+	else /* df == sd->tb */
994	+	sdt->sidef = SD_BXMIT;
995	+	sdt->st = NULL;
996	+	df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
997	+	df->comp[0].dist = sdt;
998	+	df->comp[0].func = &SDhandleTre;
999	+	/* read BSDF data */
1000	+	sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
1001	+	if (!sdata \|\| !next_token(&sdata)) {
1002	+	sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
1003	+	sd->name);
1004	+	return SDEformat;
1005	+	}
1006	+	sdt->st = load_tree_data(&sdata, ndim);
1007	+	if (sdt->st == NULL)
1008	+	return SDEformat;
1009	+	if (next_token(&sdata)) { /* check for unconsumed characters */
1010	+	sprintf(SDerrorDetail,
1011	+	"Extra characters at end of ScatteringData in '%s'",
1012	+	sd->name);
1013	+	return SDEformat;
1014	+	}
1015	+	/* flatten branches where possible */
1016	+	sdt->st = SDsimplifyTre(sdt->st);
1017	+	if (sdt->st == NULL)
1018	+	return SDEinternal;
1019	+	return get_extrema(df); /* compute global quantities */
1020	+	}
1021	+
1022	+	/* Find minimum value in tree */
1023	+	static float
1024	+	SDgetTreMin(const SDNode *st)
1025	+	{
1026	+	float vmin = FHUGE;
1027	+	int n;
1028	+
1029	+	if (st->log2GR < 0) {
1030	+	for (n = 1<<st->ndim; n--; ) {
1031	+	float v = SDgetTreMin(st->u.t[n]);
1032	+	if (v < vmin)
1033	+	vmin = v;
1034	+	}
1035	+	} else {
1036	+	for (n = 1<<(st->ndim*st->log2GR); n--; )
1037	+	if (st->u.v[n] < vmin)
1038	+	vmin = st->u.v[n];
1039	+	}
1040	+	return vmin;
1041	+	}
1042	+
1043	+	/* Subtract the given value from all tree nodes */
1044	+	static void
1045	+	SDsubtractTreVal(SDNode *st, float val)
1046	+	{
1047	+	int n;
1048	+
1049	+	if (st->log2GR < 0) {
1050	+	for (n = 1<<st->ndim; n--; )
1051	+	SDsubtractTreVal(st->u.t[n], val);
1052	+	} else {
1053	+	for (n = 1<<(st->ndim*st->log2GR); n--; )
1054	+	if ((st->u.v[n] -= val) < 0)
1055	+	st->u.v[n] = .0f;
1056	+	}
1057	+	}
1058	+
1059	+	/* Subtract minimum value from BSDF */
1060	+	static double
1061	+	subtract_min(SDNode *st)
1062	+	{
1063	+	float vmin;
1064	+	/* be sure to skip unused portion */
1065	+	if (st->ndim == 3) {
1066	+	int n;
1067	+	vmin = 1./M_PI;
1068	+	if (st->log2GR < 0) {
1069	+	for (n = 0; n < 8; n += 2) {
1070	+	float v = SDgetTreMin(st->u.t[n]);
1071	+	if (v < vmin)
1072	+	vmin = v;
1073	+	}
1074	+	} else if (st->log2GR) {
1075	+	for (n = 1 << (3*st->log2GR - 1); n--; )
1076	+	if (st->u.v[n] < vmin)
1077	+	vmin = st->u.v[n];
1078	+	} else
1079	+	vmin = st->u.v[0];
1080	+	} else /* anisotropic covers entire tree */
1081	+	vmin = SDgetTreMin(st);
1082	+
1083	+	if (vmin <= FTINY)
1084	+	return .0;
1085	+
1086	+	SDsubtractTreVal(st, vmin);
1087	+
1088	+	return M_PI * vmin; /* return hemispherical value */
1089	+	}
1090	+
1091	+	/* Extract and separate diffuse portion of BSDF */
1092	+	static void
1093	+	extract_diffuse(SDValue dv, SDSpectralDF df)
1094	+	{
1095	+	int n;
1096	+
1097	+	if (df == NULL \|\| df->ncomp <= 0) {
1098	+	dv->spec = c_dfcolor;
1099	+	dv->cieY = .0;
1100	+	return;
1101	+	}
1102	+	dv->spec = df->comp[0].cspec[0];
1103	+	dv->cieY = subtract_min(((SDTre )df->comp[0].dist).st);
1104	+	/* in case of multiple components */
1105	+	for (n = df->ncomp; --n; ) {
1106	+	double ymin = subtract_min(((SDTre )df->comp[n].dist).st);
1107	+	c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
1108	+	dv->cieY += ymin;
1109	+	}
1110	+	df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */
1111	+	/* make sure everything is set */
1112	+	c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
1113	+	}
1114	+
1115		/* Load a variable-resolution BSDF tree from an open XML file */
1116		SDError
1117		SDloadTre(SDData *sd, ezxml_t wtl)
1118		{
1119	<	return SDEsupport;
1119	>	SDError ec;
1120	>	ezxml_t wld, wdb;
1121	>	int rank;
1122	>	char *txt;
1123	>	/* basic checks and tensor rank */
1124	>	txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
1125	>	"DataDefinition"), "IncidentDataStructure"));
1126	>	if (txt == NULL \|\| !*txt) {
1127	>	sprintf(SDerrorDetail,
1128	>	"BSDF \"%s\": missing IncidentDataStructure",
1129	>	sd->name);
1130	>	return SDEformat;
1131	>	}
1132	>	if (!strcasecmp(txt, "TensorTree3"))
1133	>	rank = 3;
1134	>	else if (!strcasecmp(txt, "TensorTree4"))
1135	>	rank = 4;
1136	>	else {
1137	>	sprintf(SDerrorDetail,
1138	>	"BSDF \"%s\": unsupported IncidentDataStructure",
1139	>	sd->name);
1140	>	return SDEsupport;
1141	>	}
1142	>	/* load BSDF components */
1143	>	for (wld = ezxml_child(wtl, "WavelengthData");
1144	>	wld != NULL; wld = wld->next) {
1145	>	if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
1146	>	"Visible"))
1147	>	continue; /* just visible for now */
1148	>	for (wdb = ezxml_child(wld, "WavelengthDataBlock");
1149	>	wdb != NULL; wdb = wdb->next)
1150	>	if ((ec = load_bsdf_data(sd, wdb, rank)) != SDEnone)
1151	>	return ec;
1152	>	}
1153	>	/* separate diffuse components */
1154	>	extract_diffuse(&sd->rLambFront, sd->rf);
1155	>	extract_diffuse(&sd->rLambBack, sd->rb);
1156	>	extract_diffuse(&sd->tLamb, (sd->tf != NULL) ? sd->tf : sd->tb);
1157	>	/* return success */
1158	>	return SDEnone;
1159		}
1160
1161		/* Variable resolution BSDF methods */

Diff Legend

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

Comparing ray/src/common/bsdf_t.c (file contents): Revision 3.6 by greg, Sun Apr 24 19:39:21 2011 UTC vs. Revision 3.24 by greg, Sun Sep 2 15:33:15 2012 UTC

Diff Legend

Comparing ray/src/common/bsdf_t.c (file contents):
Revision 3.6 by greg, Sun Apr 24 19:39:21 2011 UTC vs.
Revision 3.24 by greg, Sun Sep 2 15:33:15 2012 UTC