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.6 by greg, Fri Sep 3 23:53:50 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       7               /* 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 custom 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 || !*abname)
212	+	return;
213	+	for (i = nabases; i--; )
214	+	if (!strcmp(abname, abase_list[i].name))
215	+	return;         /* 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 double
243	+	to_meters(              /* return factor to convert given unit to meters */
244	+	const char *unit
245	+	)
246	+	{
247	+	if (unit == NULL) return(1.);           /* safe assumption? */
248	+	if (!strcasecmp(unit, "Meter")) return(1.);
249	+	if (!strcasecmp(unit, "Foot")) return(.3048);
250	+	if (!strcasecmp(unit, "Inch")) return(.0254);
251	+	if (!strcasecmp(unit, "Centimeter")) return(.01);
252	+	sprintf(errmsg, "unknown dimensional unit '%s'", unit);
253	+	error(USER, errmsg);
254	+	}
255	+
256	+
257	+	static void
258	+	load_geometry(          /* load geometric dimensions and description (if any) */
259	+	struct BSDF_data *dp,
260	+	ezxml_t wdb
261	+	)
262	+	{
263	+	ezxml_t         geom;
264	+	double          cfact;
265	+	const char      *fmt, *mgfstr;
266	+
267	+	dp->dim[0] = dp->dim[1] = dp->dim[2] = 0;
268	+	dp->mgf = NULL;
269	+	if ((geom = ezxml_child(wdb, "Width")) != NULL)
270	+	dp->dim[0] = atof(ezxml_txt(geom)) *
271	+	to_meters(ezxml_attr(geom, "unit"));
272	+	if ((geom = ezxml_child(wdb, "Height")) != NULL)
273	+	dp->dim[1] = atof(ezxml_txt(geom)) *
274	+	to_meters(ezxml_attr(geom, "unit"));
275	+	if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
276	+	dp->dim[2] = atof(ezxml_txt(geom)) *
277	+	to_meters(ezxml_attr(geom, "unit"));
278	+	if ((geom = ezxml_child(wdb, "Geometry")) == NULL ||
279	+	(mgfstr = ezxml_txt(geom)) == NULL)
280	+	return;
281	+	if ((fmt = ezxml_attr(geom, "format")) != NULL &&
282	+	strcasecmp(fmt, "MGF")) {
283	+	sprintf(errmsg, "unrecognized geometry format '%s'", fmt);
284	+	error(WARNING, errmsg);
285	+	return;
286	+	}
287	+	cfact = to_meters(ezxml_attr(geom, "unit"));
288	+	dp->mgf = (char *)malloc(strlen(mgfstr)+32);
289	+	if (dp->mgf == NULL)
290	+	error(SYSTEM, "out of memory in load_geometry");
291	+	if (cfact < 0.99 || cfact > 1.01)
292	+	sprintf(dp->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
293	+	else
294	+	strcpy(dp->mgf, mgfstr);
295	+	}
296	+
297	+
298	+	static void
299		load_bsdf_data(         /* load BSDF distribution for this wavelength */
300		struct BSDF_data *dp,
301		ezxml_t wdb
#	Line 184 | Line 306 | load_bsdf_data(                /* load BSDF distribution for this wa
306		char  *sdata;
307		int  i;
308
309	<	if ((cbasis == NULL) | (rbasis == NULL)) {
309	>	if ((!cbasis || !*cbasis) | (!rbasis || !*rbasis)) {
310		error(WARNING, "missing column/row basis for BSDF");
311		return;
312		}
191	–	/* XXX need to add routines for loading in foreign bases */
313		for (i = nabases; i--; )
314		if (!strcmp(cbasis, abase_list[i].name)) {
315		dp->ninc = abase_list[i].nangles;
#	Line 199 | Line 320 | load_bsdf_data(                /* load BSDF distribution for this wa
320		break;
321		}
322		if (i < 0) {
323	<	sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
323	>	sprintf(errmsg, "undefined ColumnAngleBasis '%s'", cbasis);
324		error(WARNING, errmsg);
325		return;
326		}
#	Line 213 | Line 334 | load_bsdf_data(                /* load BSDF distribution for this wa
334		break;
335		}
336		if (i < 0) {
337	<	sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
337	>	sprintf(errmsg, "undefined RowAngleBasis '%s'", cbasis);
338		error(WARNING, errmsg);
339		return;
340		}
341		/* read BSDF data */
342		sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
343	<	if (sdata == NULL) {
343	>	if (!sdata || !*sdata) {
344		error(WARNING, "missing BSDF ScatteringData");
345		return;
346		}
#	Line 243 | Line 364 | load_bsdf_data(                /* load BSDF distribution for this wa
364		sdata++;
365		if (*sdata) {
366		sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
367	<	strlen(sdata));
367	>	(int)strlen(sdata));
368		error(WARNING, errmsg);
369		}
370		}
#	Line 254 | Line 375 | check_bsdf_data(       /* check that BSDF data is sane */
375		struct BSDF_data *dp
376		)
377		{
378	<	double *        omega_arr;
379	<	double          dom, hemi_total;
378	>	double          *omega_iarr, *omega_oarr;
379	>	double          dom, contrib, hemi_total;
380		int             nneg;
381	+	FVECT           v;
382		int             i, o;
383
384		if (dp == NULL || dp->bsdf == NULL)
385		return(0);
386	<	omega_arr = (double *)calloc(dp->nout, sizeof(double));
387	<	if (omega_arr == NULL)
386	>	omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
387	>	omega_oarr = (double *)calloc(dp->nout, sizeof(double));
388	>	if ((omega_iarr == NULL) | (omega_oarr == NULL))
389		error(SYSTEM, "out of memory in check_bsdf_data");
390	+	/* incoming projected solid angles */
391		hemi_total = .0;
392	+	for (i = dp->ninc; i--; ) {
393	+	dom = getBSDF_incohm(dp,i);
394	+	if (dom <= .0) {
395	+	error(WARNING, "zero/negative incoming solid angle");
396	+	continue;
397	+	}
398	+	if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
399	+	error(WARNING, "illegal incoming BSDF direction");
400	+	free(omega_iarr); free(omega_oarr);
401	+	return(0);
402	+	}
403	+	hemi_total += omega_iarr[i] = dom * -v[2];
404	+	}
405	+	if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
406	+	sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
407	+	100.*(hemi_total/PI - 1.));
408	+	error(WARNING, errmsg);
409	+	}
410	+	dom = PI / hemi_total;          /* fix normalization */
411	+	for (i = dp->ninc; i--; )
412	+	omega_iarr[i] *= dom;
413	+	/* outgoing projected solid angles */
414	+	hemi_total = .0;
415		for (o = dp->nout; o--; ) {
269	–	FVECT   v;
416		dom = getBSDF_outohm(dp,o);
417		if (dom <= .0) {
418	<	error(WARNING, "zero/negative solid angle");
418	>	error(WARNING, "zero/negative outgoing solid angle");
419		continue;
420		}
421		if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
422		error(WARNING, "illegal outgoing BSDF direction");
423	<	free(omega_arr);
423	>	free(omega_iarr); free(omega_oarr);
424		return(0);
425		}
426	<	hemi_total += omega_arr[o] = dom*v[2];
426	>	hemi_total += omega_oarr[o] = dom * v[2];
427		}
428		if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
429		sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
430		100.*(hemi_total/PI - 1.));
431		error(WARNING, errmsg);
432		}
433	<	dom = PI / hemi_total;          /* normalize solid angles */
433	>	dom = PI / hemi_total;          /* fix normalization */
434		for (o = dp->nout; o--; )
435	<	omega_arr[o] *= dom;
436	<	nneg = 0;
437	<	for (i = dp->ninc; i--; ) {
435	>	omega_oarr[o] *= dom;
436	>	nneg = 0;                       /* check outgoing totals */
437	>	for (i = 0; i < dp->ninc; i++) {
438		hemi_total = .0;
439		for (o = dp->nout; o--; ) {
440		double  f = BSDF_value(dp,i,o);
441	<	if (f > .0)
442	<	hemi_total += f*omega_arr[o];
443	<	else if (f < -FTINY)
444	<	++nneg;
441	>	if (f >= .0)
442	>	hemi_total += f*omega_oarr[o];
443	>	else {
444	>	nneg += (f < -FTINY);
445	>	BSDF_value(dp,i,o) = .0f;
446	>	}
447		}
448		if (hemi_total > 1.02) {
449	<	sprintf(errmsg, "BSDF direction passes %.1f%% of light",
450	<	100.*hemi_total);
449	>	sprintf(errmsg,
450	>	"incoming BSDF direction %d passes %.1f%% of light",
451	>	i, 100.*hemi_total);
452		error(WARNING, errmsg);
453		}
454		}
455	<	free(omega_arr);
456	<	if (nneg > 0) {
308	<	sprintf(errmsg, "%d negative BSDF values", nneg);
455	>	if (nneg) {
456	>	sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
457		error(WARNING, errmsg);
310	–	return(0);
458		}
459	+	/* reverse roles and check again */
460	+	for (o = 0; o < dp->nout; o++) {
461	+	hemi_total = .0;
462	+	for (i = dp->ninc; i--; )
463	+	hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
464	+
465	+	if (hemi_total > 1.02) {
466	+	sprintf(errmsg,
467	+	"outgoing BSDF direction %d collects %.1f%% of light",
468	+	o, 100.*hemi_total);
469	+	error(WARNING, errmsg);
470	+	}
471	+	}
472	+	free(omega_iarr); free(omega_oarr);
473		return(1);
474		}
475
476	+
477		struct BSDF_data *
478		load_BSDF(              /* load BSDF data from file */
479		char *fname
#	Line 348 | Line 510 | load_BSDF(             /* load BSDF data from file */
510		return(NULL);
511		}
512		wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
513	+	load_angle_basis(ezxml_child(ezxml_child(wtl,
514	+	"DataDefinition"), "AngleBasis"));
515		dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
516	+	load_geometry(dp, ezxml_child(wtl, "Material"));
517		for (wld = ezxml_child(wtl, "WavelengthData");
518		wld != NULL; wld = wld->next) {
519		if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
#	Line 379 | Line 544 | free_BSDF(             /* free BSDF data structure */
544		{
545		if (b == NULL)
546		return;
547	+	if (b->mgf != NULL)
548	+	free(b->mgf);
549		if (b->bsdf != NULL)
550		free(b->bsdf);
551		free(b);
#	Line 435 | Line 602 | r_BSDF_outvec(         /* compute random output vector at giv
602		}
603
604
438	–	#define  FEQ(a,b)       ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
439	–
605		static int
606		addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
607		char *xfarg[],
#	Line 487 | Line 652 | int
652		getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
653		MAT4 xm,
654		FVECT nrm,
655	<	UpDir ud
655	>	UpDir ud,
656	>	char *xfbuf
657		)
658		{
659		char    *xfargs[7];
660		XF      myxf;
661		FVECT   updir, xdest, ydest;
662	+	int     i;
663
664		updir[0] = updir[1] = updir[2] = 0.;
665		switch (ud) {
#	Line 523 | Line 690 | getBSDF_xfm(           /* compute BSDF orient. -> world orient.
690		fcross(ydest, nrm, xdest);
691		xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
692		copymat4(xm, myxf.xfm);
693	+	if (xfbuf == NULL)
694	+	return(1);
695	+	/* return xf arguments as well */
696	+	for (i = 0; xfargs[i] != NULL; i++) {
697	+	*xfbuf++ = ' ';
698	+	strcpy(xfbuf, xfargs[i]);
699	+	while (*xfbuf) ++xfbuf;
700	+	}
701		return(1);
702		}

Diff Legend

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