ViewVC Help

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

root/Development/ray/src/common/bsdf.c

Comparing ray/src/common/bsdf.c (file contents):
Revision 2.1 by greg, Wed Jun 17 20:41:47 2009 UTC vs.
Revision 2.3 by greg, Sat Jul 3 02:50:22 2010 UTC

#	Line 23 | Line 23 | typedef struct {
23		}       lat[MAXLATS+1];         /* latitudes */
24		} ANGLE_BASIS;
25
26	<	#define MAXABASES       3               /* limit on defined bases */
26	>	#define MAXABASES       5               /* limit on defined bases */
27
28		static ANGLE_BASIS      abase_list[MAXABASES] = {
29		{
#	Line 61 | Line 61 | static ANGLE_BASIS     abase_list[MAXABASES] = {
61
62		static int      nabases = 3;    /* current number of defined bases */
63
64	+	#define  FEQ(a,b)       ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
65
66	+	// returns the name of the given tag
67	+	#ifdef ezxml_name
68	+	#undef ezxml_name
69	+	static char *
70	+	ezxml_name(ezxml_t xml)
71	+	{
72	+	if (xml == NULL)
73	+	return(NULL);
74	+	return(xml->name);
75	+	}
76	+	#endif
77	+
78	+	// returns the given tag's character content or empty string if none
79	+	#ifdef ezxml_txt
80	+	#undef ezxml_txt
81	+	static char *
82	+	ezxml_txt(ezxml_t xml)
83	+	{
84	+	if (xml == NULL)
85	+	return("");
86	+	return(xml->txt);
87	+	}
88	+	#endif
89	+
90	+
91		static int
92		ab_getvec(              /* get vector for this angle basis index */
93		FVECT v,
#	Line 71 | Line 97 | ab_getvec(             /* get vector for this angle basis index *
97		{
98		ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
99		int     li;
100	<	double  alt, azi, d;
100	>	double  pol, azi, d;
101
102		if ((ndx < 0) | (ndx >= ab->nangles))
103		return(0);
104		for (li = 0; ndx >= ab->lat[li].nphis; li++)
105		ndx -= ab->lat[li].nphis;
106	<	alt = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
106	>	pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
107		azi = 2.*PI*ndx/ab->lat[li].nphis;
108	<	v[2] = d = cos(alt);
109	<	d = sqrt(1. - d*d);     /* sin(alt) */
108	>	v[2] = d = cos(pol);
109	>	d = sqrt(1. - d*d);     /* sin(pol) */
110		v[0] = cos(azi)*d;
111		v[1] = sin(azi)*d;
112		return(1);
#	Line 95 | Line 121 | ab_getndx(             /* get index corresponding to the given ve
121		{
122		ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
123		int     li, ndx;
124	<	double  alt, azi, d;
124	>	double  pol, azi, d;
125
126		if ((v[2] < -1.0) | (v[2] > 1.0))
127		return(-1);
128	<	alt = 180.0/PI*acos(v[2]);
128	>	pol = 180.0/PI*acos(v[2]);
129		azi = 180.0/PI*atan2(v[1], v[0]);
130		if (azi < 0.0) azi += 360.0;
131	<	for (li = 1; ab->lat[li].tmin <= alt; li++)
131	>	for (li = 1; ab->lat[li].tmin <= pol; li++)
132		if (!ab->lat[li].nphis)
133		return(-1);
134		--li;
#	Line 174 | Line 200 | ab_getndxR(            /* get index corresponding to the reverse
200
201
202		static void
203	+	load_angle_basis(       /* load BSDF angle basis */
204	+	ezxml_t wab
205	+	)
206	+	{
207	+	char    *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName"));
208	+	ezxml_t wbb;
209	+	int     i;
210	+
211	+	if (abname == NULL || !*abname)
212	+	return;
213	+	for (i = nabases; i--; )
214	+	if (!strcmp(abname, abase_list[i].name))
215	+	return;         /* XXX assume it's the same */
216	+	if (nabases >= MAXABASES)
217	+	error(INTERNAL, "too many angle bases");
218	+	strcpy(abase_list[nabases].name, abname);
219	+	abase_list[nabases].nangles = 0;
220	+	for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
221	+	wbb != NULL; i++, wbb = wbb->next) {
222	+	if (i >= MAXLATS)
223	+	error(INTERNAL, "too many latitudes in custom basis");
224	+	abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
225	+	ezxml_child(ezxml_child(wbb,
226	+	"ThetaBounds"), "UpperTheta")));
227	+	if (!i)
228	+	abase_list[nabases].lat[i].tmin =
229	+	-abase_list[nabases].lat[i+1].tmin;
230	+	else if (!FEQ(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
231	+	"ThetaBounds"), "LowerTheta"))),
232	+	abase_list[nabases].lat[i].tmin))
233	+	error(WARNING, "theta values disagree in custom basis");
234	+	abase_list[nabases].nangles +=
235	+	abase_list[nabases].lat[i].nphis =
236	+	atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
237	+	}
238	+	abase_list[nabases++].lat[i].nphis = 0;
239	+	}
240	+
241	+
242	+	static void
243		load_bsdf_data(         /* load BSDF distribution for this wavelength */
244		struct BSDF_data *dp,
245		ezxml_t wdb
#	Line 184 | Line 250 | load_bsdf_data(                /* load BSDF distribution for this wa
250		char  *sdata;
251		int  i;
252
253	<	if ((cbasis == NULL) | (rbasis == NULL)) {
253	>	if ((cbasis == NULL || !*cbasis) | (rbasis == NULL || !*rbasis)) {
254		error(WARNING, "missing column/row basis for BSDF");
255		return;
256		}
#	Line 219 | Line 285 | load_bsdf_data(                /* load BSDF distribution for this wa
285		}
286		/* read BSDF data */
287		sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
288	<	if (sdata == NULL) {
288	>	if (sdata == NULL || !*sdata) {
289		error(WARNING, "missing BSDF ScatteringData");
290		return;
291		}
#	Line 254 | Line 320 | check_bsdf_data(       /* check that BSDF data is sane */
320		struct BSDF_data *dp
321		)
322		{
323	<	double *        omega_arr;
324	<	double          dom, hemi_total;
323	>	double          *omega_iarr, *omega_oarr;
324	>	double          dom, contrib, hemi_total;
325		int             nneg;
326	+	FVECT           v;
327		int             i, o;
328
329		if (dp == NULL || dp->bsdf == NULL)
330		return(0);
331	<	omega_arr = (double *)calloc(dp->nout, sizeof(double));
332	<	if (omega_arr == NULL)
331	>	omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
332	>	omega_oarr = (double *)calloc(dp->nout, sizeof(double));
333	>	if ((omega_iarr == NULL) | (omega_oarr == NULL))
334		error(SYSTEM, "out of memory in check_bsdf_data");
335	+	/* incoming projected solid angles */
336		hemi_total = .0;
337	+	for (i = dp->ninc; i--; ) {
338	+	dom = getBSDF_incohm(dp,i);
339	+	if (dom <= .0) {
340	+	error(WARNING, "zero/negative incoming solid angle");
341	+	continue;
342	+	}
343	+	if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
344	+	error(WARNING, "illegal incoming BSDF direction");
345	+	free(omega_iarr); free(omega_oarr);
346	+	return(0);
347	+	}
348	+	hemi_total += omega_iarr[i] = dom * -v[2];
349	+	}
350	+	if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
351	+	sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
352	+	100.*(hemi_total/PI - 1.));
353	+	error(WARNING, errmsg);
354	+	}
355	+	dom = PI / hemi_total;          /* fix normalization */
356	+	for (i = dp->ninc; i--; )
357	+	omega_iarr[i] *= dom;
358	+	/* outgoing projected solid angles */
359	+	hemi_total = .0;
360		for (o = dp->nout; o--; ) {
269	–	FVECT   v;
361		dom = getBSDF_outohm(dp,o);
362		if (dom <= .0) {
363	<	error(WARNING, "zero/negative solid angle");
363	>	error(WARNING, "zero/negative outgoing solid angle");
364		continue;
365		}
366		if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
367		error(WARNING, "illegal outgoing BSDF direction");
368	<	free(omega_arr);
368	>	free(omega_iarr); free(omega_oarr);
369		return(0);
370		}
371	<	hemi_total += omega_arr[o] = dom*v[2];
371	>	hemi_total += omega_oarr[o] = dom * v[2];
372		}
373		if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
374		sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
375		100.*(hemi_total/PI - 1.));
376		error(WARNING, errmsg);
377		}
378	<	dom = PI / hemi_total;          /* normalize solid angles */
378	>	dom = PI / hemi_total;          /* fix normalization */
379		for (o = dp->nout; o--; )
380	<	omega_arr[o] *= dom;
381	<	nneg = 0;
382	<	for (i = dp->ninc; i--; ) {
380	>	omega_oarr[o] *= dom;
381	>	nneg = 0;                       /* check outgoing totals */
382	>	for (i = 0; i < dp->ninc; i++) {
383		hemi_total = .0;
384		for (o = dp->nout; o--; ) {
385		double  f = BSDF_value(dp,i,o);
386	<	if (f > .0)
387	<	hemi_total += f*omega_arr[o];
388	<	else if (f < -FTINY)
389	<	++nneg;
386	>	if (f >= .0)
387	>	hemi_total += f*omega_oarr[o];
388	>	else {
389	>	nneg += (f < -FTINY);
390	>	BSDF_value(dp,i,o) = .0f;
391	>	}
392		}
393		if (hemi_total > 1.02) {
394	<	sprintf(errmsg, "BSDF direction passes %.1f%% of light",
395	<	100.*hemi_total);
394	>	sprintf(errmsg,
395	>	"incoming BSDF direction %d passes %.1f%% of light",
396	>	i, 100.*hemi_total);
397		error(WARNING, errmsg);
398		}
399		}
400	<	free(omega_arr);
401	<	if (nneg > 0) {
308	<	sprintf(errmsg, "%d negative BSDF values", nneg);
400	>	if (nneg) {
401	>	sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
402		error(WARNING, errmsg);
310	–	return(0);
403		}
404	+	/* reverse roles and check again */
405	+	for (o = 0; o < dp->nout; o++) {
406	+	hemi_total = .0;
407	+	for (i = dp->ninc; i--; )
408	+	hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
409	+
410	+	if (hemi_total > 1.02) {
411	+	sprintf(errmsg,
412	+	"outgoing BSDF direction %d collects %.1f%% of light",
413	+	o, 100.*hemi_total);
414	+	error(WARNING, errmsg);
415	+	}
416	+	}
417	+	free(omega_iarr); free(omega_oarr);
418		return(1);
419		}
420
#	Line 348 | Line 454 | load_BSDF(             /* load BSDF data from file */
454		return(NULL);
455		}
456		wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
457	+	load_angle_basis(ezxml_child(ezxml_child(wtl,
458	+	"DataDefinition"), "AngleBasis"));
459		dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
460		for (wld = ezxml_child(wtl, "WavelengthData");
461		wld != NULL; wld = wld->next) {
#	Line 434 | Line 542 | r_BSDF_outvec(         /* compute random output vector at giv
542		return(normalize(v) != 0.0);
543		}
544
437	–
438	–	#define  FEQ(a,b)       ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
545
546		static int
547		addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)