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