[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.27 by greg, Tue Apr 23 17:25:23 2013 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 374 \| Line 499 \| *build_scaffold(float val, const double cmin, double c**
499		{
500		SDdistScaffold sp = (SDdistScaffold )cptr;
501		int wid = csiz*(double)iwmax + .5;
502	+	double revcmin[2];
503		bitmask_t bmin[2], bmax[2];
504
505	<	cmin += sp->nic; /* skip to output coords */
505	>	if (sp->rev) { /* need to reverse sense? */
506	>	revcmin[0] = 1. - cmin[0] - csiz;
507	>	revcmin[1] = 1. - cmin[1] - csiz;
508	>	cmin = revcmin;
509	>	} else {
510	>	cmin += sp->nic; /* else skip to output coords */
511	>	}
512		if (wid < sp->wmin) /* new minimum width? */
513		sp->wmin = wid;
514		if (wid > sp->wmax) /* new maximum? */
515		sp->wmax = wid;
516		if (sp->alen >= sp->nall) { /* need more space? */
517		struct outdir_s *ndarr;
518	<	sp->nall += 8192;
518	>	sp->nall += 1024;
519		ndarr = (struct outdir_s *)realloc(sp->darr,
520		sizeof(struct outdir_s)*sp->nall);
521	<	if (ndarr == NULL)
521	>	if (ndarr == NULL) {
522	>	sprintf(SDerrorDetail,
523	>	"Cannot grow scaffold to %u entries", sp->nall);
524		return -1; /* abort build */
525	+	}
526		sp->darr = ndarr;
527		}
528		/* find Hilbert entry index */
529		bmin[0] = cmin[0]*(double)iwmax + .5;
530		bmin[1] = cmin[1]*(double)iwmax + .5;
531	<	bmax[0] = bmin[0] + wid;
532	<	bmax[1] = bmin[1] + wid;
533	<	hilbert_box_vtx(2, sizeof(bitmask_t), sizeof(unsigned)*4,
534	<	1, bmin, bmax);
400	<	sp->darr[sp->alen].hent = hilbert_c2i(2, sizeof(unsigned)*4, bmin);
531	>	bmax[0] = bmin[0] + wid-1;
532	>	bmax[1] = bmin[1] + wid-1;
533	>	hilbert_box_vtx(2, sizeof(bitmask_t), iwbits, 1, bmin, bmax);
534	>	sp->darr[sp->alen].hent = hilbert_c2i(2, iwbits, bmin);
535		sp->darr[sp->alen].wid = wid;
536		sp->darr[sp->alen].bsdf = val;
537		sp->alen++; /* on to the next entry */
#	Line 408 \| Line 542 \| *build_scaffold(float val, const double cmin, double c**
542		static int
543		sscmp(const void p1, const void p2)
544		{
545	<	return (int)(((const struct outdir_s )p1).hent -
546	<	((const struct outdir_s )p2).hent);
545	>	unsigned h1 = ((const struct outdir_s )p1).hent;
546	>	unsigned h2 = ((const struct outdir_s )p2).hent;
547	>
548	>	if (h1 > h2)
549	>	return 1;
550	>	if (h1 < h2)
551	>	return -1;
552	>	return 0;
553		}
554
555		/* Create a new cumulative distribution for the given input direction */
556		static SDTreCDst *
557	<	make_cdist(const SDTre sdt, const double pos)
557	>	make_cdist(const SDTre sdt, const double invec, int rev)
558		{
419	–	const unsigned cumlmax = ~0;
559		SDdistScaffold myScaffold;
560	+	double pos[4];
561	+	int cmask;
562		SDTreCDst *cd;
563		struct outdir_s *sp;
564		double scale, cursum;
#	Line 426 \| Line 567 \| *make_cdist(const SDTre sdt, const double pos)*
567		myScaffold.wmin = iwmax;
568		myScaffold.wmax = 0;
569		myScaffold.nic = sdt->st->ndim - 2;
570	+	myScaffold.rev = rev;
571		myScaffold.alen = 0;
572	<	myScaffold.nall = 8192;
572	>	myScaffold.nall = 512;
573		myScaffold.darr = (struct outdir_s )malloc(sizeof(struct outdir_s)
574		myScaffold.nall);
575		if (myScaffold.darr == NULL)
576		return NULL;
577	+	/* set up traversal */
578	+	cmask = (1<<myScaffold.nic) - 1;
579	+	for (i = myScaffold.nic; i--; )
580	+	pos[i+2*rev] = invec[i];
581	+	cmask <<= 2*rev;
582		/* grow the distribution */
583	<	if (SDtraverseTre(sdt->st, pos, (1<<myScaffold.nic)-1,
583	>	if (SDtraverseTre(sdt->st, pos, cmask,
584		&build_scaffold, &myScaffold) < 0) {
585		free(myScaffold.darr);
586		return NULL;
#	Line 442 \| Line 589 \| *make_cdist(const SDTre sdt, const double pos)*
589		cd = (SDTreCDst *)malloc(sizeof(SDTreCDst) +
590		sizeof(cd->carr[0])*myScaffold.alen);
591		if (cd == NULL) {
592	+	sprintf(SDerrorDetail,
593	+	"Cannot allocate %u entry cumulative distribution",
594	+	myScaffold.alen);
595		free(myScaffold.darr);
596		return NULL;
597		}
598	+	cd->isodist = (myScaffold.nic == 1);
599		/* sort the distribution */
600		qsort(myScaffold.darr, cd->calen = myScaffold.alen,
601		sizeof(struct outdir_s), &sscmp);
602
603		/* record input range */
604	<	scale = (double)myScaffold.wmin / iwmax;
604	>	scale = myScaffold.wmin / (double)iwmax;
605		for (i = myScaffold.nic; i--; ) {
606	<	cd->clim[i][0] = floor(pos[i]/scale + .5) * scale;
606	>	cd->clim[i][0] = floor(pos[i+2rev]/scale) scale;
607		cd->clim[i][1] = cd->clim[i][0] + scale;
608		}
609	+	if (cd->isodist) { /* avoid issue in SDqueryTreProjSA() */
610	+	cd->clim[1][0] = cd->clim[0][0];
611	+	cd->clim[1][1] = cd->clim[0][1];
612	+	}
613		cd->max_psa = myScaffold.wmax / (double)iwmax;
614		cd->max_psa = cd->max_psa M_PI;
615	<	cd->isxmit = sdt->isxmit;
615	>	if (rev)
616	>	cd->sidef = (sdt->sidef==SD_BXMIT) ? SD_FXMIT : SD_BXMIT;
617	>	else
618	>	cd->sidef = sdt->sidef;
619		cd->cTotal = 1e-20; /* compute directional total */
620		sp = myScaffold.darr;
621		for (i = myScaffold.alen; i--; sp++)
#	Line 466 \| Line 624 \| *make_cdist(const SDTre sdt, const double pos)*
624		scale = (double)cumlmax / cd->cTotal;
625		sp = myScaffold.darr;
626		for (i = 0; i < cd->calen; i++, sp++) {
627	+	cd->carr[i].hndx = sp->hent;
628		cd->carr[i].cuml = scale*cursum + .5;
629		cursum += sp->bsdf * (double)sp->wid * sp->wid;
630		}
#	Line 483 \| Line 642 \| *SDgetTreCDist(const FVECT inVec, SDComponent sdc)**
642		{
643		const SDTre *sdt;
644		double inCoord[2];
486	–	int vflags;
645		int i;
646	+	int mode;
647		SDTreCDst cd, cdlast;
648		/* check arguments */
649		if ((inVec == NULL) \| (sdc == NULL) \|\|
650		(sdt = (SDTre *)sdc->dist) == NULL)
651		return NULL;
652	<	if (sdt->st->ndim == 3) /* isotropic BSDF? */
652	>	switch (mode = sdt->sidef) { /* check direction */
653	>	case SD_FREFL:
654	>	if (inVec[2] < 0)
655	>	return NULL;
656	>	break;
657	>	case SD_BREFL:
658	>	if (inVec[2] > 0)
659	>	return NULL;
660	>	break;
661	>	case SD_FXMIT:
662	>	if (inVec[2] < 0)
663	>	mode = SD_BXMIT;
664	>	break;
665	>	case SD_BXMIT:
666	>	if (inVec[2] > 0)
667	>	mode = SD_FXMIT;
668	>	break;
669	>	default:
670	>	return NULL;
671	>	}
672	>	if (sdt->st->ndim == 3) { /* isotropic BSDF? */
673	>	if (mode != sdt->sidef) /* XXX unhandled reciprocity */
674	>	return &SDemptyCD;
675		inCoord[0] = .5 - .5sqrt(inVec[0]inVec[0] + inVec[1]*inVec[1]);
676	<	else if (sdt->st->ndim == 4)
677	<	SDdisk2square(inCoord, -inVec[0], -inVec[1]);
678	<	else
676	>	} else if (sdt->st->ndim == 4) {
677	>	if (mode != sdt->sidef) /* use reciprocity? */
678	>	SDdisk2square(inCoord, inVec[0], inVec[1]);
679	>	else
680	>	SDdisk2square(inCoord, -inVec[0], -inVec[1]);
681	>	} else
682		return NULL; /* should be internal error */
683	+	/* quantize to avoid f.p. errors */
684	+	for (i = sdt->st->ndim - 2; i--; )
685	+	inCoord[i] = floor(inCoord[i]/quantum)quantum + .5quantum;
686		cdlast = NULL; /* check for direction in cache list */
687		for (cd = (SDTreCDst *)sdc->cdList; cd != NULL;
688	<	cdlast = cd, cd = (SDTreCDst *)cd->next) {
688	>	cdlast = cd, cd = cd->next) {
689	>	if (cd->sidef != mode)
690	>	continue;
691		for (i = sdt->st->ndim - 2; i--; )
692		if ((cd->clim[i][0] > inCoord[i]) \|
693		(inCoord[i] >= cd->clim[i][1]))
#	Line 507 \| Line 696 \| *SDgetTreCDist(const FVECT inVec, SDComponent sdc)**
696		break; /* means we have a match */
697		}
698		if (cd == NULL) /* need to create new entry? */
699	<	cdlast = cd = make_cdist(sdt, inCoord);
699	>	cdlast = cd = make_cdist(sdt, inCoord, mode != sdt->sidef);
700		if (cdlast != NULL) { /* move entry to head of cache list */
701		cdlast->next = cd->next;
702	<	cd->next = sdc->cdList;
702	>	cd->next = (SDTreCDst *)sdc->cdList;
703		sdc->cdList = (SDCDst *)cd;
704		}
705		return (SDCDst )cd; / ready to go */
#	Line 531 \| Line 720 \| *SDqueryTreProjSA(double psa, const FVECT v1, const RR**
720		const SDTre sdt = (SDTre )sdc->dist;
721		double hcube[SD_MAXDIM];
722		if (SDqueryTre(sdt, v1, v2, hcube) < 0) {
723	<	if (qflags == SDqueryVal)
724	<	*psa = M_PI;
536	<	return SDEnone;
723	>	strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA");
724	>	return SDEinternal;
725		}
726		myPSA[0] = hcube[sdt->st->ndim];
727		myPSA[1] = myPSA[0] = myPSA[0] M_PI;
#	Line 571 \| Line 759 \| SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
759		{
760		const unsigned nBitsC = 4*sizeof(bitmask_t);
761		const unsigned nExtraBits = 8*(sizeof(bitmask_t)-sizeof(unsigned));
574	–	const unsigned maxval = ~0;
762		const SDTreCDst cd = (const SDTreCDst )cdp;
763	<	const unsigned target = randX*maxval;
763	>	const unsigned target = randX*cumlmax;
764		bitmask_t hndx, hcoord[2];
765	<	double gpos[3];
765	>	double gpos[3], rotangle;
766		int i, iupper, ilower;
767		/* check arguments */
768		if ((ioVec == NULL) \| (cd == NULL))
769		return SDEargument;
770	+	if (!cd->sidef)
771	+	return SDEnone; /* XXX should never happen */
772	+	if (ioVec[2] > 0) {
773	+	if ((cd->sidef != SD_FREFL) & (cd->sidef != SD_FXMIT))
774	+	return SDEargument;
775	+	} else if ((cd->sidef != SD_BREFL) & (cd->sidef != SD_BXMIT))
776	+	return SDEargument;
777		/* binary search to find position */
778		ilower = 0; iupper = cd->calen;
779		while ((i = (iupper + ilower) >> 1) != ilower)
780	<	if ((long)target >= (long)cd->carr[i].cuml)
780	>	if (target >= cd->carr[i].cuml)
781		ilower = i;
782		else
783		iupper = i;
784		/* localize random position */
785	<	randX = (randX*maxval - cd->carr[ilower].cuml) /
785	>	randX = (randX*cumlmax - cd->carr[ilower].cuml) /
786		(double)(cd->carr[iupper].cuml - cd->carr[ilower].cuml);
787		/* index in longer Hilbert curve */
788		hndx = (randXcd->carr[iupper].hndx + (1.-randX)cd->carr[ilower].hndx)
#	Line 599 \| Line 793 \| SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
793		gpos[i] = ((double)hcoord[i] + rand()*(1./(RAND_MAX+.5))) /
794		(double)((bitmask_t)1 << nBitsC);
795		SDsquare2disk(gpos, gpos[0], gpos[1]);
796	+	/* compute Z-coordinate */
797		gpos[2] = 1. - gpos[0]gpos[0] - gpos[1]gpos[1];
798		if (gpos[2] > 0) /* paranoia, I hope */
799		gpos[2] = sqrt(gpos[2]);
800	<	if (ioVec[2] > 0 ^ !cd->isxmit)
800	>	/* emit from back? */
801	>	if ((cd->sidef == SD_BREFL) \| (cd->sidef == SD_FXMIT))
802		gpos[2] = -gpos[2];
803	<	VCOPY(ioVec, gpos);
803	>	if (cd->isodist) { /* rotate isotropic result */
804	>	rotangle = atan2(-ioVec[1],-ioVec[0]);
805	>	VCOPY(ioVec, gpos);
806	>	spinvector(ioVec, ioVec, zvec, rotangle);
807	>	} else
808	>	VCOPY(ioVec, gpos);
809		return SDEnone;
810		}
811
812	+	/* Advance pointer to the next non-white character in the string (or nul) */
813	+	static int
814	+	next_token(char **spp)
815	+	{
816	+	while (isspace(**spp))
817	+	++*spp;
818	+	return **spp;
819	+	}
820	+
821	+	/* Advance pointer past matching token (or any token if c==0) */
822	+	#define eat_token(spp,c) (next_token(spp)==(c) ^ !(c) ? ((spp))++ : 0)
823	+
824	+	/* Count words from this point in string to '}' */
825	+	static int
826	+	count_values(char *cp)
827	+	{
828	+	int n = 0;
829	+
830	+	while (next_token(&cp) != '}' && *cp) {
831	+	while (!isspace(cp) & (cp != ',') & (*cp != '}'))
832	+	if (!*++cp)
833	+	break;
834	+	++n;
835	+	eat_token(&cp, ',');
836	+	}
837	+	return n;
838	+	}
839	+
840	+	/* Load an array of real numbers, returning total */
841	+	static int
842	+	load_values(char *spp, float va, int n)
843	+	{
844	+	float *v = va;
845	+	char *svnext;
846	+
847	+	while (n-- > 0 && (svnext = fskip(*spp)) != NULL) {
848	+	v++ = atof(spp);
849	+	*spp = svnext;
850	+	eat_token(spp, ',');
851	+	}
852	+	return v - va;
853	+	}
854	+
855	+	/* Load BSDF tree data */
856	+	static SDNode *
857	+	load_tree_data(char **spp, int nd)
858	+	{
859	+	SDNode *st;
860	+	int n;
861	+
862	+	if (!eat_token(spp, '{')) {
863	+	strcpy(SDerrorDetail, "Missing '{' in tensor tree");
864	+	return NULL;
865	+	}
866	+	if (next_token(spp) == '{') { /* tree branches */
867	+	st = SDnewNode(nd, -1);
868	+	if (st == NULL)
869	+	return NULL;
870	+	for (n = 0; n < 1<<nd; n++)
871	+	if ((st->u.t[n] = load_tree_data(spp, nd)) == NULL) {
872	+	SDfreeTre(st);
873	+	return NULL;
874	+	}
875	+	} else { /* else load value grid */
876	+	int bsiz;
877	+	n = count_values(spp); / see how big the grid is */
878	+	for (bsiz = 0; bsiz < 8*sizeof(size_t); bsiz += nd)
879	+	if (1<<bsiz == n)
880	+	break;
881	+	if (bsiz >= 8*sizeof(size_t)) {
882	+	strcpy(SDerrorDetail, "Illegal value count in tensor tree");
883	+	return NULL;
884	+	}
885	+	st = SDnewNode(nd, bsiz/nd);
886	+	if (st == NULL)
887	+	return NULL;
888	+	if (load_values(spp, st->u.v, n) != n) {
889	+	strcpy(SDerrorDetail, "Real format error in tensor tree");
890	+	SDfreeTre(st);
891	+	return NULL;
892	+	}
893	+	}
894	+	if (!eat_token(spp, '}')) {
895	+	strcpy(SDerrorDetail, "Missing '}' in tensor tree");
896	+	SDfreeTre(st);
897	+	return NULL;
898	+	}
899	+	eat_token(spp, ',');
900	+	return st;
901	+	}
902	+
903	+	/* Compute min. proj. solid angle and max. direct hemispherical scattering */
904	+	static SDError
905	+	get_extrema(SDSpectralDF *df)
906	+	{
907	+	SDNode st = ((SDTre *)df->comp[0].dist).st;
908	+	double stepWidth, dhemi, bmin[4], bmax[4];
909	+
910	+	stepWidth = SDsmallestLeaf(st);
911	+	if (quantum > stepWidth) /* adjust quantization factor */
912	+	quantum = stepWidth;
913	+	df->minProjSA = M_PIstepWidthstepWidth;
914	+	if (stepWidth < .03125)
915	+	stepWidth = .03125; /* 1/32 resolution good enough */
916	+	df->maxHemi = .0;
917	+	if (st->ndim == 3) { /* isotropic BSDF */
918	+	bmin[1] = bmin[2] = .0;
919	+	bmax[1] = bmax[2] = 1.;
920	+	for (bmin[0] = .0; bmin[0] < .5-FTINY; bmin[0] += stepWidth) {
921	+	bmax[0] = bmin[0] + stepWidth;
922	+	dhemi = SDavgTreBox(st, bmin, bmax);
923	+	if (dhemi > df->maxHemi)
924	+	df->maxHemi = dhemi;
925	+	}
926	+	} else if (st->ndim == 4) { /* anisotropic BSDF */
927	+	bmin[2] = bmin[3] = .0;
928	+	bmax[2] = bmax[3] = 1.;
929	+	for (bmin[0] = .0; bmin[0] < 1.-FTINY; bmin[0] += stepWidth) {
930	+	bmax[0] = bmin[0] + stepWidth;
931	+	for (bmin[1] = .0; bmin[1] < 1.-FTINY; bmin[1] += stepWidth) {
932	+	bmax[1] = bmin[1] + stepWidth;
933	+	dhemi = SDavgTreBox(st, bmin, bmax);
934	+	if (dhemi > df->maxHemi)
935	+	df->maxHemi = dhemi;
936	+	}
937	+	}
938	+	} else
939	+	return SDEinternal;
940	+	/* correct hemispherical value */
941	+	df->maxHemi *= M_PI;
942	+	return SDEnone;
943	+	}
944	+
945	+	/* Load BSDF distribution for this wavelength */
946	+	static SDError
947	+	load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
948	+	{
949	+	SDSpectralDF *df;
950	+	SDTre *sdt;
951	+	char *sdata;
952	+	/* allocate BSDF component */
953	+	sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
954	+	if (!sdata)
955	+	return SDEnone;
956	+	/*
957	+	* Remember that front and back are reversed from WINDOW 6 orientations
958	+	*/
959	+	if (!strcasecmp(sdata, "Transmission Front")) {
960	+	if (sd->tb != NULL)
961	+	SDfreeSpectralDF(sd->tb);
962	+	if ((sd->tb = SDnewSpectralDF(1)) == NULL)
963	+	return SDEmemory;
964	+	df = sd->tb;
965	+	} else if (!strcasecmp(sdata, "Transmission Back")) {
966	+	if (sd->tf != NULL)
967	+	SDfreeSpectralDF(sd->tf);
968	+	if ((sd->tf = SDnewSpectralDF(1)) == NULL)
969	+	return SDEmemory;
970	+	df = sd->tf;
971	+	} else if (!strcasecmp(sdata, "Reflection Front")) {
972	+	if (sd->rb != NULL)
973	+	SDfreeSpectralDF(sd->rb);
974	+	if ((sd->rb = SDnewSpectralDF(1)) == NULL)
975	+	return SDEmemory;
976	+	df = sd->rb;
977	+	} else if (!strcasecmp(sdata, "Reflection Back")) {
978	+	if (sd->rf != NULL)
979	+	SDfreeSpectralDF(sd->rf);
980	+	if ((sd->rf = SDnewSpectralDF(1)) == NULL)
981	+	return SDEmemory;
982	+	df = sd->rf;
983	+	} else
984	+	return SDEnone;
985	+	/* XXX should also check "ScatteringDataType" for consistency? */
986	+	/* get angle bases */
987	+	sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis"));
988	+	if (!sdata \|\| strcasecmp(sdata, "LBNL/Shirley-Chiu")) {
989	+	sprintf(SDerrorDetail, "%s angle basis for BSDF '%s'",
990	+	!sdata ? "Missing" : "Unsupported", sd->name);
991	+	return !sdata ? SDEformat : SDEsupport;
992	+	}
993	+	/* allocate BSDF tree */
994	+	sdt = (SDTre *)malloc(sizeof(SDTre));
995	+	if (sdt == NULL)
996	+	return SDEmemory;
997	+	if (df == sd->rf)
998	+	sdt->sidef = SD_FREFL;
999	+	else if (df == sd->rb)
1000	+	sdt->sidef = SD_BREFL;
1001	+	else if (df == sd->tf)
1002	+	sdt->sidef = SD_FXMIT;
1003	+	else /* df == sd->tb */
1004	+	sdt->sidef = SD_BXMIT;
1005	+	sdt->st = NULL;
1006	+	df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
1007	+	df->comp[0].dist = sdt;
1008	+	df->comp[0].func = &SDhandleTre;
1009	+	/* read BSDF data */
1010	+	sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
1011	+	if (!sdata \|\| !next_token(&sdata)) {
1012	+	sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
1013	+	sd->name);
1014	+	return SDEformat;
1015	+	}
1016	+	sdt->st = load_tree_data(&sdata, ndim);
1017	+	if (sdt->st == NULL)
1018	+	return SDEformat;
1019	+	if (next_token(&sdata)) { /* check for unconsumed characters */
1020	+	sprintf(SDerrorDetail,
1021	+	"Extra characters at end of ScatteringData in '%s'",
1022	+	sd->name);
1023	+	return SDEformat;
1024	+	}
1025	+	/* flatten branches where possible */
1026	+	sdt->st = SDsimplifyTre(sdt->st);
1027	+	if (sdt->st == NULL)
1028	+	return SDEinternal;
1029	+	return get_extrema(df); /* compute global quantities */
1030	+	}
1031	+
1032	+	/* Find minimum value in tree */
1033	+	static float
1034	+	SDgetTreMin(const SDNode *st)
1035	+	{
1036	+	float vmin = FHUGE;
1037	+	int n;
1038	+
1039	+	if (st->log2GR < 0) {
1040	+	for (n = 1<<st->ndim; n--; ) {
1041	+	float v = SDgetTreMin(st->u.t[n]);
1042	+	if (v < vmin)
1043	+	vmin = v;
1044	+	}
1045	+	} else {
1046	+	for (n = 1<<(st->ndim*st->log2GR); n--; )
1047	+	if (st->u.v[n] < vmin)
1048	+	vmin = st->u.v[n];
1049	+	}
1050	+	return vmin;
1051	+	}
1052	+
1053	+	/* Subtract the given value from all tree nodes */
1054	+	static void
1055	+	SDsubtractTreVal(SDNode *st, float val)
1056	+	{
1057	+	int n;
1058	+
1059	+	if (st->log2GR < 0) {
1060	+	for (n = 1<<st->ndim; n--; )
1061	+	SDsubtractTreVal(st->u.t[n], val);
1062	+	} else {
1063	+	for (n = 1<<(st->ndim*st->log2GR); n--; )
1064	+	if ((st->u.v[n] -= val) < 0)
1065	+	st->u.v[n] = .0f;
1066	+	}
1067	+	}
1068	+
1069	+	/* Subtract minimum value from BSDF */
1070	+	static double
1071	+	subtract_min(SDNode *st)
1072	+	{
1073	+	float vmin;
1074	+	/* be sure to skip unused portion */
1075	+	if (st->ndim == 3) {
1076	+	int n;
1077	+	vmin = 1./M_PI;
1078	+	if (st->log2GR < 0) {
1079	+	for (n = 0; n < 8; n += 2) {
1080	+	float v = SDgetTreMin(st->u.t[n]);
1081	+	if (v < vmin)
1082	+	vmin = v;
1083	+	}
1084	+	} else if (st->log2GR) {
1085	+	for (n = 1 << (3*st->log2GR - 1); n--; )
1086	+	if (st->u.v[n] < vmin)
1087	+	vmin = st->u.v[n];
1088	+	} else
1089	+	vmin = st->u.v[0];
1090	+	} else /* anisotropic covers entire tree */
1091	+	vmin = SDgetTreMin(st);
1092	+
1093	+	if (vmin <= FTINY)
1094	+	return .0;
1095	+
1096	+	SDsubtractTreVal(st, vmin);
1097	+
1098	+	return M_PI * vmin; /* return hemispherical value */
1099	+	}
1100	+
1101	+	/* Extract and separate diffuse portion of BSDF */
1102	+	static void
1103	+	extract_diffuse(SDValue dv, SDSpectralDF df)
1104	+	{
1105	+	int n;
1106	+
1107	+	if (df == NULL \|\| df->ncomp <= 0) {
1108	+	dv->spec = c_dfcolor;
1109	+	dv->cieY = .0;
1110	+	return;
1111	+	}
1112	+	dv->spec = df->comp[0].cspec[0];
1113	+	dv->cieY = subtract_min(((SDTre )df->comp[0].dist).st);
1114	+	/* in case of multiple components */
1115	+	for (n = df->ncomp; --n; ) {
1116	+	double ymin = subtract_min(((SDTre )df->comp[n].dist).st);
1117	+	c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
1118	+	dv->cieY += ymin;
1119	+	}
1120	+	df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */
1121	+	/* make sure everything is set */
1122	+	c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
1123	+	}
1124	+
1125		/* Load a variable-resolution BSDF tree from an open XML file */
1126		SDError
1127		SDloadTre(SDData *sd, ezxml_t wtl)
1128		{
1129	<	return SDEsupport;
1129	>	SDError ec;
1130	>	ezxml_t wld, wdb;
1131	>	int rank;
1132	>	char *txt;
1133	>	/* basic checks and tensor rank */
1134	>	txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
1135	>	"DataDefinition"), "IncidentDataStructure"));
1136	>	if (txt == NULL \|\| !*txt) {
1137	>	sprintf(SDerrorDetail,
1138	>	"BSDF \"%s\": missing IncidentDataStructure",
1139	>	sd->name);
1140	>	return SDEformat;
1141	>	}
1142	>	if (!strcasecmp(txt, "TensorTree3"))
1143	>	rank = 3;
1144	>	else if (!strcasecmp(txt, "TensorTree4"))
1145	>	rank = 4;
1146	>	else {
1147	>	sprintf(SDerrorDetail,
1148	>	"BSDF \"%s\": unsupported IncidentDataStructure",
1149	>	sd->name);
1150	>	return SDEsupport;
1151	>	}
1152	>	/* load BSDF components */
1153	>	for (wld = ezxml_child(wtl, "WavelengthData");
1154	>	wld != NULL; wld = wld->next) {
1155	>	if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
1156	>	"Visible"))
1157	>	continue; /* just visible for now */
1158	>	for (wdb = ezxml_child(wld, "WavelengthDataBlock");
1159	>	wdb != NULL; wdb = wdb->next)
1160	>	if ((ec = load_bsdf_data(sd, wdb, rank)) != SDEnone)
1161	>	return ec;
1162	>	}
1163	>	/* separate diffuse components */
1164	>	extract_diffuse(&sd->rLambFront, sd->rf);
1165	>	extract_diffuse(&sd->rLambBack, sd->rb);
1166	>	if (sd->tf != NULL)
1167	>	extract_diffuse(&sd->tLamb, sd->tf);
1168	>	if (sd->tb != NULL)
1169	>	extract_diffuse(&sd->tLamb, sd->tb);
1170	>	/* return success */
1171	>	return SDEnone;
1172		}
1173
1174		/* 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.27 by greg, Tue Apr 23 17:25:23 2013 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.27 by greg, Tue Apr 23 17:25:23 2013 UTC