ViewVC Help

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

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

Comparing ray/src/common/bsdf.c (file contents):
Revision 2.9 by greg, Mon Jan 3 19:27:00 2011 UTC vs.
Revision 2.43 by greg, Sat Oct 13 20:15:43 2012 UTC

#	Line 2 | Line 2
2		static const char RCSid[] = "$Id$";
3		#endif
4		/*
5	<	* Routines for handling BSDF data
5	>	*  bsdf.c
6	>	*
7	>	*  Definitions for bidirectional scattering distribution functions.
8	>	*
9	>	*  Created by Greg Ward on 1/10/11.
10	>	*
11		*/
12
13	<	#include "standard.h"
14	<	#include "bsdf.h"
15	<	#include "paths.h"
16	<	#include "ezxml.h"
13	>	#define _USE_MATH_DEFINES
14	>	#include <stdio.h>
15	>	#include <stdlib.h>
16	>	#include <string.h>
17	>	#include <math.h>
18		#include <ctype.h>
19	+	#include "ezxml.h"
20	+	#include "hilbert.h"
21	+	#include "bsdf.h"
22	+	#include "bsdf_m.h"
23	+	#include "bsdf_t.h"
24
25	<	#define MAXLATS         46              /* maximum number of latitudes */
25	>	/* English ASCII strings corresponding to ennumerated errors */
26	>	const char              *SDerrorEnglish[] = {
27	>	"No error",
28	>	"Memory error",
29	>	"File input/output error",
30	>	"File format error",
31	>	"Illegal argument",
32	>	"Invalid data",
33	>	"Unsupported feature",
34	>	"Internal program error",
35	>	"Unknown error"
36	>	};
37
38	<	/* BSDF angle specification */
39	<	typedef struct {
18	<	char    name[64];               /* basis name */
19	<	int     nangles;                /* total number of directions */
20	<	struct {
21	<	float   tmin;                   /* starting theta */
22	<	short   nphis;                  /* number of phis (0 term) */
23	<	}       lat[MAXLATS+1];         /* latitudes */
24	<	} ANGLE_BASIS;
38	>	/* Additional information on last error (ASCII English) */
39	>	char                    SDerrorDetail[256];
40
41	<	#define MAXABASES       7               /* limit on defined bases */
41	>	/* Empty distribution for getCDist() calls that fail for some reason */
42	>	const SDCDst            SDemptyCD;
43
44	<	static ANGLE_BASIS      abase_list[MAXABASES] = {
45	<	{
30	<	"LBNL/Klems Full", 145,
31	<	{ {-5., 1},
32	<	{5., 8},
33	<	{15., 16},
34	<	{25., 20},
35	<	{35., 24},
36	<	{45., 24},
37	<	{55., 24},
38	<	{65., 16},
39	<	{75., 12},
40	<	{90., 0} }
41	<	}, {
42	<	"LBNL/Klems Half", 73,
43	<	{ {-6.5, 1},
44	<	{6.5, 8},
45	<	{19.5, 12},
46	<	{32.5, 16},
47	<	{46.5, 20},
48	<	{61.5, 12},
49	<	{76.5, 4},
50	<	{90., 0} }
51	<	}, {
52	<	"LBNL/Klems Quarter", 41,
53	<	{ {-9., 1},
54	<	{9., 8},
55	<	{27., 12},
56	<	{46., 12},
57	<	{66., 8},
58	<	{90., 0} }
59	<	}
60	<	};
44	>	/* Cache of loaded BSDFs */
45	>	struct SDCache_s        *SDcacheList = NULL;
46
47	<	static int      nabases = 3;    /* current number of defined bases */
47	>	/* Retain BSDFs in cache list */
48	>	int                     SDretainSet = SDretainNone;
49
50	<	#define  FEQ(a,b)       ((a)-(b) <= 1e-6 && (b)-(a) <= 1e-6)
51	<
52	<	static int
67	<	fequal(double a, double b)
50	>	/* Report any error to the indicated stream (in English) */
51	>	SDError
52	>	SDreportEnglish(SDError ec, FILE *fp)
53		{
54	<	if (b != .0)
55	<	a = a/b - 1.;
56	<	return((a <= 1e-6) & (a >= -1e-6));
54	>	if (!ec)
55	>	return SDEnone;
56	>	if ((ec < SDEnone) | (ec > SDEunknown)) {
57	>	SDerrorDetail[0] = '\0';
58	>	ec = SDEunknown;
59	>	}
60	>	if (fp == NULL)
61	>	return ec;
62	>	fputs(SDerrorEnglish[ec], fp);
63	>	if (SDerrorDetail[0]) {
64	>	fputs(": ", fp);
65	>	fputs(SDerrorDetail, fp);
66	>	}
67	>	fputc('\n', fp);
68	>	if (fp != stderr)
69	>	fflush(fp);
70	>	return ec;
71		}
72
73	<	// returns the name of the given tag
74	<	#ifdef ezxml_name
75	<	#undef ezxml_name
76	<	static char *
78	<	ezxml_name(ezxml_t xml)
73	>	static double
74	>	to_meters(              /* return factor to convert given unit to meters */
75	>	const char *unit
76	>	)
77		{
78	<	if (xml == NULL)
79	<	return(NULL);
80	<	return(xml->name);
78	>	if (unit == NULL) return(1.);           /* safe assumption? */
79	>	if (!strcasecmp(unit, "Meter")) return(1.);
80	>	if (!strcasecmp(unit, "Foot")) return(.3048);
81	>	if (!strcasecmp(unit, "Inch")) return(.0254);
82	>	if (!strcasecmp(unit, "Centimeter")) return(.01);
83	>	if (!strcasecmp(unit, "Millimeter")) return(.001);
84	>	sprintf(SDerrorDetail, "Unknown dimensional unit '%s'", unit);
85	>	return(-1.);
86		}
84	–	#endif
87
88	<	// returns the given tag's character content or empty string if none
89	<	#ifdef ezxml_txt
90	<	#undef ezxml_txt
89	<	static char *
90	<	ezxml_txt(ezxml_t xml)
88	>	/* Load geometric dimensions and description (if any) */
89	>	static SDError
90	>	SDloadGeometry(SDData *sd, ezxml_t wtl)
91		{
92	<	if (xml == NULL)
93	<	return("");
94	<	return(xml->txt);
95	<	}
96	<	#endif
92	>	ezxml_t         node, matl, geom;
93	>	double          cfact;
94	>	const char      *fmt = NULL, *mgfstr;
95
96	<
97	<	static int
98	<	ab_getvec(              /* get vector for this angle basis index */
99	<	FVECT v,
100	<	int ndx,
101	<	void *p
102	<	)
103	<	{
104	<	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
105	<	int     li;
106	<	double  pol, azi, d;
107	<
108	<	if ((ndx < 0) | (ndx >= ab->nangles))
109	<	return(0);
110	<	for (li = 0; ndx >= ab->lat[li].nphis; li++)
111	<	ndx -= ab->lat[li].nphis;
112	<	pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
113	<	azi = 2.*PI*ndx/ab->lat[li].nphis;
114	<	v[2] = d = cos(pol);
115	<	d = sqrt(1. - d*d);     /* sin(pol) */
116	<	v[0] = cos(azi)*d;
117	<	v[1] = sin(azi)*d;
118	<	return(1);
96	>	SDerrorDetail[0] = '\0';
97	>	sd->matn[0] = '\0'; sd->makr[0] = '\0';
98	>	sd->dim[0] = sd->dim[1] = sd->dim[2] = 0;
99	>	matl = ezxml_child(wtl, "Material");
100	>	if (matl != NULL) {                     /* get material info. */
101	>	if ((node = ezxml_child(matl, "Name")) != NULL) {
102	>	strncpy(sd->matn, ezxml_txt(node), SDnameLn);
103	>	if (sd->matn[SDnameLn-1])
104	>	strcpy(sd->matn+(SDnameLn-4), "...");
105	>	}
106	>	if ((node = ezxml_child(matl, "Manufacturer")) != NULL) {
107	>	strncpy(sd->makr, ezxml_txt(node), SDnameLn);
108	>	if (sd->makr[SDnameLn-1])
109	>	strcpy(sd->makr+(SDnameLn-4), "...");
110	>	}
111	>	if ((node = ezxml_child(matl, "Width")) != NULL)
112	>	sd->dim[0] = atof(ezxml_txt(node)) *
113	>	to_meters(ezxml_attr(node, "unit"));
114	>	if ((node = ezxml_child(matl, "Height")) != NULL)
115	>	sd->dim[1] = atof(ezxml_txt(node)) *
116	>	to_meters(ezxml_attr(node, "unit"));
117	>	if ((node = ezxml_child(matl, "Thickness")) != NULL)
118	>	sd->dim[2] = atof(ezxml_txt(node)) *
119	>	to_meters(ezxml_attr(node, "unit"));
120	>	if ((sd->dim[0] < 0) | (sd->dim[1] < 0) | (sd->dim[2] < 0)) {
121	>	if (!SDerrorDetail[0])
122	>	sprintf(SDerrorDetail, "Negative dimension in \"%s\"",
123	>	sd->name);
124	>	return SDEdata;
125	>	}
126	>	}
127	>	sd->mgf = NULL;
128	>	geom = ezxml_child(wtl, "Geometry");
129	>	if (geom == NULL)                       /* no actual geometry? */
130	>	return SDEnone;
131	>	fmt = ezxml_attr(geom, "format");
132	>	if (fmt != NULL && strcasecmp(fmt, "MGF")) {
133	>	sprintf(SDerrorDetail,
134	>	"Unrecognized geometry format '%s' in \"%s\"",
135	>	fmt, sd->name);
136	>	return SDEsupport;
137	>	}
138	>	if ((node = ezxml_child(geom, "MGFblock")) == NULL ||
139	>	(mgfstr = ezxml_txt(node)) == NULL)
140	>	return SDEnone;
141	>	while (isspace(*mgfstr))
142	>	++mgfstr;
143	>	if (!*mgfstr)
144	>	return SDEnone;
145	>	cfact = to_meters(ezxml_attr(node, "unit"));
146	>	if (cfact <= 0)
147	>	return SDEformat;
148	>	sd->mgf = (char *)malloc(strlen(mgfstr)+32);
149	>	if (sd->mgf == NULL) {
150	>	strcpy(SDerrorDetail, "Out of memory in SDloadGeometry");
151	>	return SDEmemory;
152	>	}
153	>	if (cfact < 0.99 || cfact > 1.01)
154	>	sprintf(sd->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
155	>	else
156	>	strcpy(sd->mgf, mgfstr);
157	>	return SDEnone;
158		}
159
160	<
161	<	static int
162	<	ab_getndx(              /* get index corresponding to the given vector */
126	<	FVECT v,
127	<	void *p
128	<	)
160	>	/* Load a BSDF struct from the given file (free first and keep name) */
161	>	SDError
162	>	SDloadFile(SDData *sd, const char *fname)
163		{
164	<	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
165	<	int     li, ndx;
132	<	double  pol, azi, d;
164	>	SDError         lastErr;
165	>	ezxml_t         fl, wtl;
166
167	<	if ((v[2] < -1.0) | (v[2] > 1.0))
168	<	return(-1);
169	<	pol = 180.0/PI*acos(v[2]);
170	<	azi = 180.0/PI*atan2(v[1], v[0]);
171	<	if (azi < 0.0) azi += 360.0;
172	<	for (li = 1; ab->lat[li].tmin <= pol; li++)
173	<	if (!ab->lat[li].nphis)
174	<	return(-1);
175	<	--li;
176	<	ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
177	<	if (ndx >= ab->lat[li].nphis) ndx = 0;
178	<	while (li--)
179	<	ndx += ab->lat[li].nphis;
180	<	return(ndx);
167	>	if ((sd == NULL) | (fname == NULL || !*fname))
168	>	return SDEargument;
169	>	/* free old data, keeping name */
170	>	SDfreeBSDF(sd);
171	>	/* parse XML file */
172	>	fl = ezxml_parse_file(fname);
173	>	if (fl == NULL) {
174	>	sprintf(SDerrorDetail, "Cannot open BSDF \"%s\"", fname);
175	>	return SDEfile;
176	>	}
177	>	if (ezxml_error(fl)[0]) {
178	>	sprintf(SDerrorDetail, "BSDF \"%s\" %s", fname, ezxml_error(fl));
179	>	ezxml_free(fl);
180	>	return SDEformat;
181	>	}
182	>	if (strcmp(ezxml_name(fl), "WindowElement")) {
183	>	sprintf(SDerrorDetail,
184	>	"BSDF \"%s\": top level node not 'WindowElement'",
185	>	sd->name);
186	>	ezxml_free(fl);
187	>	return SDEformat;
188	>	}
189	>	wtl = ezxml_child(fl, "FileType");
190	>	if (wtl != NULL && strcmp(ezxml_txt(wtl), "BSDF")) {
191	>	sprintf(SDerrorDetail,
192	>	"XML \"%s\": wrong FileType (must be 'BSDF')",
193	>	sd->name);
194	>	ezxml_free(fl);
195	>	return SDEformat;
196	>	}
197	>	wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
198	>	if (wtl == NULL) {
199	>	sprintf(SDerrorDetail, "BSDF \"%s\": no optical layers'",
200	>	sd->name);
201	>	ezxml_free(fl);
202	>	return SDEformat;
203	>	}
204	>	/* load geometry if present */
205	>	lastErr = SDloadGeometry(sd, wtl);
206	>	if (lastErr) {
207	>	ezxml_free(fl);
208	>	return lastErr;
209	>	}
210	>	/* try loading variable resolution data */
211	>	lastErr = SDloadTre(sd, wtl);
212	>	/* check our result */
213	>	if (lastErr == SDEsupport)      /* try matrix BSDF if not tree data */
214	>	lastErr = SDloadMtx(sd, wtl);
215	>
216	>	/* done with XML file */
217	>	ezxml_free(fl);
218	>
219	>	if (lastErr) {          /* was there a load error? */
220	>	SDfreeBSDF(sd);
221	>	return lastErr;
222	>	}
223	>	/* remove any insignificant components */
224	>	if (sd->rf != NULL && sd->rf->maxHemi <= .001) {
225	>	SDfreeSpectralDF(sd->rf); sd->rf = NULL;
226	>	}
227	>	if (sd->rb != NULL && sd->rb->maxHemi <= .001) {
228	>	SDfreeSpectralDF(sd->rb); sd->rb = NULL;
229	>	}
230	>	if (sd->tf != NULL && sd->tf->maxHemi <= .001) {
231	>	SDfreeSpectralDF(sd->tf); sd->tf = NULL;
232	>	}
233	>	if (sd->tb != NULL && sd->tb->maxHemi <= .001) {
234	>	SDfreeSpectralDF(sd->tb); sd->tb = NULL;
235	>	}
236	>	/* return success */
237	>	return SDEnone;
238		}
239
240	<
241	<	static double
242	<	ab_getohm(              /* get solid angle for this angle basis index */
153	<	int ndx,
154	<	void *p
155	<	)
240	>	/* Allocate new spectral distribution function */
241	>	SDSpectralDF *
242	>	SDnewSpectralDF(int nc)
243		{
244	<	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
158	<	int     li;
159	<	double  theta, theta1;
244	>	SDSpectralDF    *df;
245
246	<	if ((ndx < 0) | (ndx >= ab->nangles))
247	<	return(0);
248	<	for (li = 0; ndx >= ab->lat[li].nphis; li++)
164	<	ndx -= ab->lat[li].nphis;
165	<	theta1 = PI/180. * ab->lat[li+1].tmin;
166	<	if (ab->lat[li].nphis == 1) {           /* special case */
167	<	if (ab->lat[li].tmin > FTINY)
168	<	error(USER, "unsupported BSDF coordinate system");
169	<	return(2.*PI*(1. - cos(theta1)));
246	>	if (nc <= 0) {
247	>	strcpy(SDerrorDetail, "Zero component spectral DF request");
248	>	return NULL;
249		}
250	<	theta = PI/180. * ab->lat[li].tmin;
251	<	return(2.*PI*(cos(theta) - cos(theta1))/(double)ab->lat[li].nphis);
250	>	df = (SDSpectralDF *)malloc(sizeof(SDSpectralDF) +
251	>	(nc-1)*sizeof(SDComponent));
252	>	if (df == NULL) {
253	>	sprintf(SDerrorDetail,
254	>	"Cannot allocate %d component spectral DF", nc);
255	>	return NULL;
256	>	}
257	>	df->minProjSA = .0;
258	>	df->maxHemi = .0;
259	>	df->ncomp = nc;
260	>	memset(df->comp, 0, nc*sizeof(SDComponent));
261	>	return df;
262		}
263
264	<
265	<	static int
266	<	ab_getvecR(             /* get reverse vector for this angle basis index */
178	<	FVECT v,
179	<	int ndx,
180	<	void *p
181	<	)
264	>	/* Add component(s) to spectral distribution function */
265	>	SDSpectralDF *
266	>	SDaddComponent(SDSpectralDF *odf, int nadd)
267		{
268	<	if (!ab_getvec(v, ndx, p))
184	<	return(0);
268	>	SDSpectralDF    *df;
269
270	<	v[0] = -v[0];
271	<	v[1] = -v[1];
272	<	v[2] = -v[2];
273	<
274	<	return(1);
270	>	if (odf == NULL)
271	>	return SDnewSpectralDF(nadd);
272	>	if (nadd <= 0)
273	>	return odf;
274	>	df = (SDSpectralDF *)realloc(odf, sizeof(SDSpectralDF) +
275	>	(odf->ncomp+nadd-1)*sizeof(SDComponent));
276	>	if (df == NULL) {
277	>	sprintf(SDerrorDetail,
278	>	"Cannot add %d component(s) to spectral DF", nadd);
279	>	SDfreeSpectralDF(odf);
280	>	return NULL;
281	>	}
282	>	memset(df->comp+df->ncomp, 0, nadd*sizeof(SDComponent));
283	>	df->ncomp += nadd;
284	>	return df;
285		}
286
287	<
288	<	static int
289	<	ab_getndxR(             /* get index corresponding to the reverse vector */
196	<	FVECT v,
197	<	void *p
198	<	)
287	>	/* Free cached cumulative distributions for BSDF component */
288	>	void
289	>	SDfreeCumulativeCache(SDSpectralDF *df)
290		{
291	<	FVECT  v2;
292	<
202	<	v2[0] = -v[0];
203	<	v2[1] = -v[1];
204	<	v2[2] = -v[2];
291	>	int     n;
292	>	SDCDst  *cdp;
293
294	<	return ab_getndx(v2, p);
294	>	if (df == NULL)
295	>	return;
296	>	for (n = df->ncomp; n-- > 0; )
297	>	while ((cdp = df->comp[n].cdList) != NULL) {
298	>	df->comp[n].cdList = cdp->next;
299	>	free(cdp);
300	>	}
301		}
302
303	<
304	<	static void
305	<	load_angle_basis(       /* load custom BSDF angle basis */
212	<	ezxml_t wab
213	<	)
303	>	/* Free a spectral distribution function */
304	>	void
305	>	SDfreeSpectralDF(SDSpectralDF *df)
306		{
307	<	char    *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName"));
308	<	ezxml_t wbb;
309	<	int     i;
218	<
219	<	if (!abname || !*abname)
307	>	int     n;
308	>
309	>	if (df == NULL)
310		return;
311	<	for (i = nabases; i--; )
312	<	if (!strcmp(abname, abase_list[i].name))
313	<	return;         /* assume it's the same */
314	<	if (nabases >= MAXABASES)
315	<	error(INTERNAL, "too many angle bases");
226	<	strcpy(abase_list[nabases].name, abname);
227	<	abase_list[nabases].nangles = 0;
228	<	for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
229	<	wbb != NULL; i++, wbb = wbb->next) {
230	<	if (i >= MAXLATS)
231	<	error(INTERNAL, "too many latitudes in custom basis");
232	<	abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
233	<	ezxml_child(ezxml_child(wbb,
234	<	"ThetaBounds"), "UpperTheta")));
235	<	if (!i)
236	<	abase_list[nabases].lat[i].tmin =
237	<	-abase_list[nabases].lat[i+1].tmin;
238	<	else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
239	<	"ThetaBounds"), "LowerTheta"))),
240	<	abase_list[nabases].lat[i].tmin))
241	<	error(WARNING, "theta values disagree in custom basis");
242	<	abase_list[nabases].nangles +=
243	<	abase_list[nabases].lat[i].nphis =
244	<	atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
245	<	}
246	<	abase_list[nabases++].lat[i].nphis = 0;
311	>	SDfreeCumulativeCache(df);
312	>	for (n = df->ncomp; n-- > 0; )
313	>	if (df->comp[n].dist != NULL)
314	>	(*df->comp[n].func->freeSC)(df->comp[n].dist);
315	>	free(df);
316		}
317
318	<
319	<	static double
320	<	to_meters(              /* return factor to convert given unit to meters */
252	<	const char *unit
253	<	)
318	>	/* Shorten file path to useable BSDF name, removing suffix */
319	>	void
320	>	SDclipName(char *res, const char *fname)
321		{
322	<	if (unit == NULL) return(1.);           /* safe assumption? */
323	<	if (!strcasecmp(unit, "Meter")) return(1.);
324	<	if (!strcasecmp(unit, "Foot")) return(.3048);
325	<	if (!strcasecmp(unit, "Inch")) return(.0254);
326	<	if (!strcasecmp(unit, "Centimeter")) return(.01);
327	<	if (!strcasecmp(unit, "Millimeter")) return(.001);
328	<	sprintf(errmsg, "unknown dimensional unit '%s'", unit);
329	<	error(USER, errmsg);
322	>	const char      *cp, *dot = NULL;
323	>
324	>	for (cp = fname; *cp; cp++)
325	>	if (*cp == '.')
326	>	dot = cp;
327	>	if ((dot == NULL) | (dot < fname+2))
328	>	dot = cp;
329	>	if (dot - fname >= SDnameLn)
330	>	fname = dot - SDnameLn + 1;
331	>	while (fname < dot)
332	>	*res++ = *fname++;
333	>	*res = '\0';
334		}
335
336	<
337	<	static void
338	<	load_geometry(          /* load geometric dimensions and description (if any) */
268	<	struct BSDF_data *dp,
269	<	ezxml_t wdb
270	<	)
336	>	/* Initialize an unused BSDF struct (simply clears to zeroes) */
337	>	void
338	>	SDclearBSDF(SDData *sd, const char *fname)
339		{
340	<	ezxml_t         geom;
273	<	double          cfact;
274	<	const char      *fmt, *mgfstr;
275	<
276	<	dp->dim[0] = dp->dim[1] = dp->dim[2] = 0;
277	<	dp->mgf = NULL;
278	<	if ((geom = ezxml_child(wdb, "Width")) != NULL)
279	<	dp->dim[0] = atof(ezxml_txt(geom)) *
280	<	to_meters(ezxml_attr(geom, "unit"));
281	<	if ((geom = ezxml_child(wdb, "Height")) != NULL)
282	<	dp->dim[1] = atof(ezxml_txt(geom)) *
283	<	to_meters(ezxml_attr(geom, "unit"));
284	<	if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
285	<	dp->dim[2] = atof(ezxml_txt(geom)) *
286	<	to_meters(ezxml_attr(geom, "unit"));
287	<	if ((geom = ezxml_child(wdb, "Geometry")) == NULL ||
288	<	(mgfstr = ezxml_txt(geom)) == NULL)
340	>	if (sd == NULL)
341		return;
342	<	if ((fmt = ezxml_attr(geom, "format")) != NULL &&
343	<	strcasecmp(fmt, "MGF")) {
292	<	sprintf(errmsg, "unrecognized geometry format '%s'", fmt);
293	<	error(WARNING, errmsg);
342	>	memset(sd, 0, sizeof(SDData));
343	>	if (fname == NULL)
344		return;
345	<	}
296	<	cfact = to_meters(ezxml_attr(geom, "unit"));
297	<	dp->mgf = (char *)malloc(strlen(mgfstr)+32);
298	<	if (dp->mgf == NULL)
299	<	error(SYSTEM, "out of memory in load_geometry");
300	<	if (cfact < 0.99 || cfact > 1.01)
301	<	sprintf(dp->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
302	<	else
303	<	strcpy(dp->mgf, mgfstr);
345	>	SDclipName(sd->name, fname);
346		}
347
348	<
349	<	static void
350	<	load_bsdf_data(         /* load BSDF distribution for this wavelength */
309	<	struct BSDF_data *dp,
310	<	ezxml_t wdb
311	<	)
348	>	/* Free data associated with BSDF struct */
349	>	void
350	>	SDfreeBSDF(SDData *sd)
351		{
352	<	char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
314	<	char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
315	<	char  *sdata;
316	<	int  i;
317	<
318	<	if ((!cbasis || !*cbasis) | (!rbasis || !*rbasis)) {
319	<	error(WARNING, "missing column/row basis for BSDF");
352	>	if (sd == NULL)
353		return;
354	+	if (sd->mgf != NULL) {
355	+	free(sd->mgf);
356	+	sd->mgf = NULL;
357		}
358	<	for (i = nabases; i--; )
359	<	if (!strcmp(cbasis, abase_list[i].name)) {
360	<	dp->ninc = abase_list[i].nangles;
325	<	dp->ib_priv = (void *)&abase_list[i];
326	<	dp->ib_vec = ab_getvecR;
327	<	dp->ib_ndx = ab_getndxR;
328	<	dp->ib_ohm = ab_getohm;
329	<	break;
330	<	}
331	<	if (i < 0) {
332	<	sprintf(errmsg, "undefined ColumnAngleBasis '%s'", cbasis);
333	<	error(WARNING, errmsg);
334	<	return;
358	>	if (sd->rf != NULL) {
359	>	SDfreeSpectralDF(sd->rf);
360	>	sd->rf = NULL;
361		}
362	<	for (i = nabases; i--; )
363	<	if (!strcmp(rbasis, abase_list[i].name)) {
364	<	dp->nout = abase_list[i].nangles;
339	<	dp->ob_priv = (void *)&abase_list[i];
340	<	dp->ob_vec = ab_getvec;
341	<	dp->ob_ndx = ab_getndx;
342	<	dp->ob_ohm = ab_getohm;
343	<	break;
344	<	}
345	<	if (i < 0) {
346	<	sprintf(errmsg, "undefined RowAngleBasis '%s'", cbasis);
347	<	error(WARNING, errmsg);
348	<	return;
362	>	if (sd->rb != NULL) {
363	>	SDfreeSpectralDF(sd->rb);
364	>	sd->rb = NULL;
365		}
366	<	/* read BSDF data */
367	<	sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
368	<	if (!sdata || !*sdata) {
353	<	error(WARNING, "missing BSDF ScatteringData");
354	<	return;
366	>	if (sd->tf != NULL) {
367	>	SDfreeSpectralDF(sd->tf);
368	>	sd->tf = NULL;
369		}
370	<	dp->bsdf = (float *)malloc(sizeof(float)*dp->ninc*dp->nout);
371	<	if (dp->bsdf == NULL)
372	<	error(SYSTEM, "out of memory in load_bsdf_data");
359	<	for (i = 0; i < dp->ninc*dp->nout; i++) {
360	<	char  *sdnext = fskip(sdata);
361	<	if (sdnext == NULL) {
362	<	error(WARNING, "bad/missing BSDF ScatteringData");
363	<	free(dp->bsdf); dp->bsdf = NULL;
364	<	return;
365	<	}
366	<	while (*sdnext && isspace(*sdnext))
367	<	sdnext++;
368	<	if (*sdnext == ',') sdnext++;
369	<	dp->bsdf[i] = atof(sdata);
370	<	sdata = sdnext;
370	>	if (sd->tb != NULL) {
371	>	SDfreeSpectralDF(sd->tb);
372	>	sd->tb = NULL;
373		}
374	<	while (isspace(*sdata))
375	<	sdata++;
376	<	if (*sdata) {
377	<	sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
378	<	(int)strlen(sdata));
379	<	error(WARNING, errmsg);
378	<	}
374	>	sd->rLambFront.cieY = .0;
375	>	sd->rLambFront.spec.flags = 0;
376	>	sd->rLambBack.cieY = .0;
377	>	sd->rLambBack.spec.flags = 0;
378	>	sd->tLamb.cieY = .0;
379	>	sd->tLamb.spec.flags = 0;
380		}
381
382	<
383	<	static int
384	<	check_bsdf_data(        /* check that BSDF data is sane */
384	<	struct BSDF_data *dp
385	<	)
382	>	/* Find writeable BSDF by name, or allocate new cache entry if absent */
383	>	SDData *
384	>	SDgetCache(const char *bname)
385		{
386	<	double          *omega_iarr, *omega_oarr;
387	<	double          dom, contrib, hemi_total, full_total;
389	<	int             nneg;
390	<	FVECT           v;
391	<	int             i, o;
386	>	struct SDCache_s        *sdl;
387	>	char                    sdnam[SDnameLn];
388
389	<	if (dp == NULL || dp->bsdf == NULL)
390	<	return(0);
391	<	omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
392	<	omega_oarr = (double *)calloc(dp->nout, sizeof(double));
393	<	if ((omega_iarr == NULL) | (omega_oarr == NULL))
394	<	error(SYSTEM, "out of memory in check_bsdf_data");
395	<	/* incoming projected solid angles */
396	<	hemi_total = .0;
401	<	for (i = dp->ninc; i--; ) {
402	<	dom = getBSDF_incohm(dp,i);
403	<	if (dom <= .0) {
404	<	error(WARNING, "zero/negative incoming solid angle");
405	<	continue;
389	>	if (bname == NULL)
390	>	return NULL;
391	>
392	>	SDclipName(sdnam, bname);
393	>	for (sdl = SDcacheList; sdl != NULL; sdl = sdl->next)
394	>	if (!strcmp(sdl->bsdf.name, sdnam)) {
395	>	sdl->refcnt++;
396	>	return &sdl->bsdf;
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];
398	>
399	>	sdl = (struct SDCache_s *)calloc(1, sizeof(struct SDCache_s));
400	>	if (sdl == NULL)
401	>	return NULL;
402	>
403	>	strcpy(sdl->bsdf.name, sdnam);
404	>	sdl->next = SDcacheList;
405	>	SDcacheList = sdl;
406	>
407	>	sdl->refcnt = 1;
408	>	return &sdl->bsdf;
409	>	}
410	>
411	>	/* Get loaded BSDF from cache (or load and cache it on first call) */
412	>	/* Report any problem to stderr and return NULL on failure */
413	>	const SDData *
414	>	SDcacheFile(const char *fname)
415	>	{
416	>	SDData          *sd;
417	>	SDError         ec;
418	>
419	>	if (fname == NULL || !*fname)
420	>	return NULL;
421	>	SDerrorDetail[0] = '\0';
422	>	if ((sd = SDgetCache(fname)) == NULL) {
423	>	SDreportEnglish(SDEmemory, stderr);
424	>	return NULL;
425		}
426	<	if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
427	<	sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
428	<	100.*(hemi_total/PI - 1.));
429	<	error(WARNING, errmsg);
426	>	if (!SDisLoaded(sd) && (ec = SDloadFile(sd, fname))) {
427	>	SDreportEnglish(ec, stderr);
428	>	SDfreeCache(sd);
429	>	return NULL;
430		}
431	<	dom = PI / hemi_total;          /* fix normalization */
432	<	for (i = dp->ninc; i--; )
421	<	omega_iarr[i] *= dom;
422	<	/* outgoing projected solid angles */
423	<	hemi_total = .0;
424	<	for (o = dp->nout; o--; ) {
425	<	dom = getBSDF_outohm(dp,o);
426	<	if (dom <= .0) {
427	<	error(WARNING, "zero/negative outgoing solid angle");
428	<	continue;
429	<	}
430	<	if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
431	<	error(WARNING, "illegal outgoing BSDF direction");
432	<	free(omega_iarr); free(omega_oarr);
433	<	return(0);
434	<	}
435	<	hemi_total += omega_oarr[o] = dom * v[2];
436	<	}
437	<	if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
438	<	sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
439	<	100.*(hemi_total/PI - 1.));
440	<	error(WARNING, errmsg);
441	<	}
442	<	dom = PI / hemi_total;          /* fix normalization */
443	<	for (o = dp->nout; o--; )
444	<	omega_oarr[o] *= dom;
445	<	nneg = 0;                       /* check outgoing totals */
446	<	for (i = 0; i < dp->ninc; i++) {
447	<	hemi_total = .0;
448	<	for (o = dp->nout; o--; ) {
449	<	double  f = BSDF_value(dp,i,o);
450	<	if (f >= .0)
451	<	hemi_total += f*omega_oarr[o];
452	<	else {
453	<	nneg += (f < -FTINY);
454	<	BSDF_value(dp,i,o) = .0f;
455	<	}
456	<	}
457	<	if (hemi_total > 1.01) {
458	<	sprintf(errmsg,
459	<	"incoming BSDF direction %d passes %.1f%% of light",
460	<	i, 100.*hemi_total);
461	<	error(WARNING, errmsg);
462	<	}
463	<	}
464	<	if (nneg) {
465	<	sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
466	<	error(WARNING, errmsg);
467	<	}
468	<	full_total = .0;                /* reverse roles and check again */
469	<	for (o = 0; o < dp->nout; o++) {
470	<	hemi_total = .0;
471	<	for (i = dp->ninc; i--; )
472	<	hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
431	>	return sd;
432	>	}
433
434	<	if (hemi_total > 1.01) {
435	<	sprintf(errmsg,
436	<	"outgoing BSDF direction %d collects %.1f%% of light",
437	<	o, 100.*hemi_total);
438	<	error(WARNING, errmsg);
434	>	/* Free a BSDF from our cache (clear all if NULL) */
435	>	void
436	>	SDfreeCache(const SDData *sd)
437	>	{
438	>	struct SDCache_s        *sdl, *sdLast = NULL;
439	>
440	>	if (sd == NULL) {               /* free entire list */
441	>	while ((sdl = SDcacheList) != NULL) {
442	>	SDcacheList = sdl->next;
443	>	SDfreeBSDF(&sdl->bsdf);
444	>	free(sdl);
445		}
446	<	full_total += hemi_total*omega_oarr[o];
446	>	return;
447		}
448	<	full_total /= PI;
449	<	if (full_total > 1.00001) {
450	<	sprintf(errmsg, "BSDF transfers %.4f%% of light",
451	<	100.*full_total);
452	<	error(WARNING, errmsg);
448	>	for (sdl = SDcacheList; sdl != NULL; sdl = (sdLast=sdl)->next)
449	>	if (&sdl->bsdf == sd)
450	>	break;
451	>	if (sdl == NULL || (sdl->refcnt -= (sdl->refcnt > 0)))
452	>	return;                 /* missing or still in use */
453	>	/* keep unreferenced data? */
454	>	if (SDisLoaded(sd) && SDretainSet) {
455	>	if (SDretainSet == SDretainAll)
456	>	return;         /* keep everything */
457	>	/* else free cumulative data */
458	>	SDfreeCumulativeCache(sd->rf);
459	>	SDfreeCumulativeCache(sd->rb);
460	>	SDfreeCumulativeCache(sd->tf);
461	>	SDfreeCumulativeCache(sd->tb);
462	>	return;
463		}
464	<	free(omega_iarr); free(omega_oarr);
465	<	return(1);
464	>	/* remove from list and free */
465	>	if (sdLast == NULL)
466	>	SDcacheList = sdl->next;
467	>	else
468	>	sdLast->next = sdl->next;
469	>	SDfreeBSDF(&sdl->bsdf);
470	>	free(sdl);
471		}
472
473	<
474	<	struct BSDF_data *
475	<	load_BSDF(              /* load BSDF data from file */
495	<	char *fname
496	<	)
473	>	/* Sample an individual BSDF component */
474	>	SDError
475	>	SDsampComponent(SDValue *sv, FVECT ioVec, double randX, SDComponent *sdc)
476		{
477	<	char                    *path;
478	<	ezxml_t                 fl, wtl, wld, wdb;
479	<	struct BSDF_data        *dp;
480	<
481	<	path = getpath(fname, getrlibpath(), R_OK);
482	<	if (path == NULL) {
483	<	sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
484	<	error(WARNING, errmsg);
485	<	return(NULL);
477	>	float           coef[SDmaxCh];
478	>	SDError         ec;
479	>	FVECT           inVec;
480	>	const SDCDst    *cd;
481	>	double          d;
482	>	int             n;
483	>	/* check arguments */
484	>	if ((sv == NULL) | (ioVec == NULL) | (sdc == NULL))
485	>	return SDEargument;
486	>	/* get cumulative distribution */
487	>	VCOPY(inVec, ioVec);
488	>	cd = (*sdc->func->getCDist)(inVec, sdc);
489	>	if (cd == NULL)
490	>	return SDEmemory;
491	>	if (cd->cTotal <= 1e-6) {       /* anything to sample? */
492	>	sv->spec = c_dfcolor;
493	>	sv->cieY = .0;
494	>	memset(ioVec, 0, 3*sizeof(double));
495	>	return SDEnone;
496		}
497	<	fl = ezxml_parse_file(path);
498	<	if (fl == NULL) {
499	<	sprintf(errmsg, "cannot open BSDF \"%s\"", path);
500	<	error(WARNING, errmsg);
501	<	return(NULL);
497	>	sv->cieY = cd->cTotal;
498	>	/* compute sample direction */
499	>	ec = (*sdc->func->sampCDist)(ioVec, randX, cd);
500	>	if (ec)
501	>	return ec;
502	>	/* get BSDF color */
503	>	n = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
504	>	if (n <= 0) {
505	>	strcpy(SDerrorDetail, "BSDF sample value error");
506	>	return SDEinternal;
507		}
508	<	if (ezxml_error(fl)[0]) {
509	<	sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
510	<	error(WARNING, errmsg);
511	<	ezxml_free(fl);
512	<	return(NULL);
508	>	sv->spec = sdc->cspec[0];
509	>	d = coef[0];
510	>	while (--n) {
511	>	c_cmix(&sv->spec, d, &sv->spec, coef[n], &sdc->cspec[n]);
512	>	d += coef[n];
513		}
514	<	if (strcmp(ezxml_name(fl), "WindowElement")) {
515	<	sprintf(errmsg,
516	<	"BSDF \"%s\": top level node not 'WindowElement'",
523	<	path);
524	<	error(WARNING, errmsg);
525	<	ezxml_free(fl);
526	<	return(NULL);
527	<	}
528	<	wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
529	<	load_angle_basis(ezxml_child(ezxml_child(wtl,
530	<	"DataDefinition"), "AngleBasis"));
531	<	dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
532	<	load_geometry(dp, ezxml_child(wtl, "Material"));
533	<	for (wld = ezxml_child(wtl, "WavelengthData");
534	<	wld != NULL; wld = wld->next) {
535	<	if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
536	<	continue;
537	<	wdb = ezxml_child(wld, "WavelengthDataBlock");
538	<	if (wdb == NULL) continue;
539	<	if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
540	<	"Transmission Front"))
541	<	continue;
542	<	load_bsdf_data(dp, wdb);        /* load front BTDF */
543	<	break;                          /* ignore the rest */
544	<	}
545	<	ezxml_free(fl);                         /* done with XML file */
546	<	if (!check_bsdf_data(dp)) {
547	<	sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
548	<	error(WARNING, errmsg);
549	<	free_BSDF(dp);
550	<	dp = NULL;
551	<	}
552	<	return(dp);
514	>	/* make sure everything is set */
515	>	c_ccvt(&sv->spec, C_CSXY+C_CSSPEC);
516	>	return SDEnone;
517		}
518
519	+	#define MS_MAXDIM       15
520
521	+	/* Convert 1-dimensional random variable to N-dimensional */
522		void
523	<	free_BSDF(              /* free BSDF data structure */
558	<	struct BSDF_data *b
559	<	)
523	>	SDmultiSamp(double t[], int n, double randX)
524		{
525	<	if (b == NULL)
525	>	unsigned        nBits;
526	>	double          scale;
527	>	bitmask_t       ndx, coord[MS_MAXDIM];
528	>
529	>	if (n <= 0)                     /* check corner cases */
530		return;
531	<	if (b->mgf != NULL)
532	<	free(b->mgf);
533	<	if (b->bsdf != NULL)
534	<	free(b->bsdf);
535	<	free(b);
531	>	if (randX < 0) randX = 0;
532	>	else if (randX >= 1.) randX = 0.999999999999999;
533	>	if (n == 1) {
534	>	t[0] = randX;
535	>	return;
536	>	}
537	>	while (n > MS_MAXDIM)           /* punt for higher dimensions */
538	>	t[--n] = rand()*(1./(RAND_MAX+.5));
539	>	nBits = (8*sizeof(bitmask_t) - 1) / n;
540	>	ndx = randX * (double)((bitmask_t)1 << (nBits*n));
541	>	/* get coordinate on Hilbert curve */
542	>	hilbert_i2c(n, nBits, ndx, coord);
543	>	/* convert back to [0,1) range */
544	>	scale = 1. / (double)((bitmask_t)1 << nBits);
545	>	while (n--)
546	>	t[n] = scale * ((double)coord[n] + rand()*(1./(RAND_MAX+.5)));
547		}
548
549	+	#undef MS_MAXDIM
550
551	<	int
552	<	r_BSDF_incvec(          /* compute random input vector at given location */
553	<	FVECT v,
574	<	struct BSDF_data *b,
575	<	int i,
576	<	double rv,
577	<	MAT4 xm
578	<	)
551	>	/* Generate diffuse hemispherical sample */
552	>	static void
553	>	SDdiffuseSamp(FVECT outVec, int outFront, double randX)
554		{
555	<	FVECT   pert;
556	<	double  rad;
557	<	int     j;
558	<
559	<	if (!getBSDF_incvec(v, b, i))
560	<	return(0);
561	<	rad = sqrt(getBSDF_incohm(b, i) / PI);
562	<	multisamp(pert, 3, rv);
588	<	for (j = 0; j < 3; j++)
589	<	v[j] += rad*(2.*pert[j] - 1.);
590	<	if (xm != NULL)
591	<	multv3(v, v, xm);
592	<	return(normalize(v) != 0.0);
555	>	/* convert to position on hemisphere */
556	>	SDmultiSamp(outVec, 2, randX);
557	>	SDsquare2disk(outVec, outVec[0], outVec[1]);
558	>	outVec[2] = 1. - outVec[0]*outVec[0] - outVec[1]*outVec[1];
559	>	if (outVec[2] > 0)              /* a bit of paranoia */
560	>	outVec[2] = sqrt(outVec[2]);
561	>	if (!outFront)                  /* going out back? */
562	>	outVec[2] = -outVec[2];
563		}
564
565	+	/* Query projected solid angle coverage for non-diffuse BSDF direction */
566	+	SDError
567	+	SDsizeBSDF(double *projSA, const FVECT v1, const RREAL *v2,
568	+	int qflags, const SDData *sd)
569	+	{
570	+	SDSpectralDF    *rdf, *tdf;
571	+	SDError         ec;
572	+	int             i;
573	+	/* check arguments */
574	+	if ((projSA == NULL) | (v1 == NULL) | (sd == NULL))
575	+	return SDEargument;
576	+	/* initialize extrema */
577	+	switch (qflags) {
578	+	case SDqueryMax:
579	+	projSA[0] = .0;
580	+	break;
581	+	case SDqueryMin+SDqueryMax:
582	+	projSA[1] = .0;
583	+	/* fall through */
584	+	case SDqueryMin:
585	+	projSA[0] = 10.;
586	+	break;
587	+	case 0:
588	+	return SDEargument;
589	+	}
590	+	if (v1[2] > 0) {                /* front surface query? */
591	+	rdf = sd->rf;
592	+	tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
593	+	} else {
594	+	rdf = sd->rb;
595	+	tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
596	+	}
597	+	if (v2 != NULL)                 /* bidirectional? */
598	+	if (v1[2] > 0 ^ v2[2] > 0)
599	+	rdf = NULL;
600	+	else
601	+	tdf = NULL;
602	+	ec = SDEdata;                   /* run through components */
603	+	for (i = (rdf==NULL) ? 0 : rdf->ncomp; i--; ) {
604	+	ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2,
605	+	qflags, &rdf->comp[i]);
606	+	if (ec)
607	+	return ec;
608	+	}
609	+	for (i = (tdf==NULL) ? 0 : tdf->ncomp; i--; ) {
610	+	ec = (*tdf->comp[i].func->queryProjSA)(projSA, v1, v2,
611	+	qflags, &tdf->comp[i]);
612	+	if (ec)
613	+	return ec;
614	+	}
615	+	if (ec) {                       /* all diffuse? */
616	+	projSA[0] = M_PI;
617	+	if (qflags == SDqueryMin+SDqueryMax)
618	+	projSA[1] = M_PI;
619	+	}
620	+	return SDEnone;
621	+	}
622
623	<	int
624	<	r_BSDF_outvec(          /* compute random output vector at given location */
625	<	FVECT v,
599	<	struct BSDF_data *b,
600	<	int o,
601	<	double rv,
602	<	MAT4 xm
603	<	)
623	>	/* Return BSDF for the given incident and scattered ray vectors */
624	>	SDError
625	>	SDevalBSDF(SDValue *sv, const FVECT outVec, const FVECT inVec, const SDData *sd)
626		{
627	<	FVECT   pert;
628	<	double  rad;
629	<	int     j;
630	<
631	<	if (!getBSDF_outvec(v, b, o))
632	<	return(0);
633	<	rad = sqrt(getBSDF_outohm(b, o) / PI);
634	<	multisamp(pert, 3, rv);
635	<	for (j = 0; j < 3; j++)
636	<	v[j] += rad*(2.*pert[j] - 1.);
637	<	if (xm != NULL)
638	<	multv3(v, v, xm);
639	<	return(normalize(v) != 0.0);
627	>	int             inFront, outFront;
628	>	SDSpectralDF    *sdf;
629	>	float           coef[SDmaxCh];
630	>	int             nch, i;
631	>	/* check arguments */
632	>	if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL))
633	>	return SDEargument;
634	>	/* whose side are we on? */
635	>	inFront = (inVec[2] > 0);
636	>	outFront = (outVec[2] > 0);
637	>	/* start with diffuse portion */
638	>	if (inFront & outFront) {
639	>	*sv = sd->rLambFront;
640	>	sdf = sd->rf;
641	>	} else if (!(inFront | outFront)) {
642	>	*sv = sd->rLambBack;
643	>	sdf = sd->rb;
644	>	} else if (inFront) {
645	>	*sv = sd->tLamb;
646	>	sdf = (sd->tf != NULL) ? sd->tf : sd->tb;
647	>	} else /* inBack */ {
648	>	*sv = sd->tLamb;
649	>	sdf = (sd->tb != NULL) ? sd->tb : sd->tf;
650	>	}
651	>	sv->cieY *= 1./M_PI;
652	>	/* add non-diffuse components */
653	>	i = (sdf != NULL) ? sdf->ncomp : 0;
654	>	while (i-- > 0) {
655	>	nch = (*sdf->comp[i].func->getBSDFs)(coef, outVec, inVec,
656	>	&sdf->comp[i]);
657	>	while (nch-- > 0) {
658	>	c_cmix(&sv->spec, sv->cieY, &sv->spec,
659	>	coef[nch], &sdf->comp[i].cspec[nch]);
660	>	sv->cieY += coef[nch];
661	>	}
662	>	}
663	>	/* make sure everything is set */
664	>	c_ccvt(&sv->spec, C_CSXY+C_CSSPEC);
665	>	return SDEnone;
666		}
667
668	<
669	<	static int
670	<	addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
623	<	char *xfarg[],
624	<	FVECT xp,
625	<	FVECT yp,
626	<	FVECT zp
627	<	)
668	>	/* Compute directional hemispherical scattering at this incident angle */
669	>	double
670	>	SDdirectHemi(const FVECT inVec, int sflags, const SDData *sd)
671		{
672	<	static char     bufs[3][16];
673	<	int     bn = 0;
674	<	char    **xfp = xfarg;
675	<	double  theta;
672	>	double          hsum;
673	>	SDSpectralDF    *rdf, *tdf;
674	>	const SDCDst    *cd;
675	>	int             i;
676	>	/* check arguments */
677	>	if ((inVec == NULL) | (sd == NULL))
678	>	return .0;
679	>	/* gather diffuse components */
680	>	if (inVec[2] > 0) {
681	>	hsum = sd->rLambFront.cieY;
682	>	rdf = sd->rf;
683	>	tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
684	>	} else /* !inFront */ {
685	>	hsum = sd->rLambBack.cieY;
686	>	rdf = sd->rb;
687	>	tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
688	>	}
689	>	if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
690	>	hsum = .0;
691	>	if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
692	>	hsum += sd->tLamb.cieY;
693	>	/* gather non-diffuse components */
694	>	i = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
695	>	(rdf != NULL)) ? rdf->ncomp : 0;
696	>	while (i-- > 0) {               /* non-diffuse reflection */
697	>	cd = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]);
698	>	if (cd != NULL)
699	>	hsum += cd->cTotal;
700	>	}
701	>	i = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
702	>	(tdf != NULL)) ? tdf->ncomp : 0;
703	>	while (i-- > 0) {               /* non-diffuse transmission */
704	>	cd = (*tdf->comp[i].func->getCDist)(inVec, &tdf->comp[i]);
705	>	if (cd != NULL)
706	>	hsum += cd->cTotal;
707	>	}
708	>	return hsum;
709	>	}
710
711	<	if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
712	<	/* Special case for X' along Z-axis */
713	<	theta = -atan2(yp[0], yp[1]);
714	<	*xfp++ = "-ry";
715	<	*xfp++ = xp[2] < 0.0 ? "90" : "-90";
716	<	*xfp++ = "-rz";
717	<	sprintf(bufs[bn], "%f", theta*(180./PI));
718	<	*xfp++ = bufs[bn++];
719	<	return(xfp - xfarg);
711	>	/* Sample BSDF direction based on the given random variable */
712	>	SDError
713	>	SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int sflags, const SDData *sd)
714	>	{
715	>	SDError         ec;
716	>	FVECT           inVec;
717	>	int             inFront;
718	>	SDSpectralDF    *rdf, *tdf;
719	>	double          rdiff;
720	>	float           coef[SDmaxCh];
721	>	int             i, j, n, nr;
722	>	SDComponent     *sdc;
723	>	const SDCDst    **cdarr = NULL;
724	>	/* check arguments */
725	>	if ((sv == NULL) | (ioVec == NULL) | (sd == NULL) |
726	>	(randX < 0) | (randX >= 1.))
727	>	return SDEargument;
728	>	/* whose side are we on? */
729	>	VCOPY(inVec, ioVec);
730	>	inFront = (inVec[2] > 0);
731	>	/* remember diffuse portions */
732	>	if (inFront) {
733	>	*sv = sd->rLambFront;
734	>	rdf = sd->rf;
735	>	tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
736	>	} else /* !inFront */ {
737	>	*sv = sd->rLambBack;
738	>	rdf = sd->rb;
739	>	tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
740		}
741	<	theta = atan2(yp[2], zp[2]);
742	<	if (!FEQ(theta,0.0)) {
743	<	*xfp++ = "-rx";
744	<	sprintf(bufs[bn], "%f", theta*(180./PI));
745	<	*xfp++ = bufs[bn++];
741	>	if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
742	>	sv->cieY = .0;
743	>	rdiff = sv->cieY;
744	>	if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
745	>	sv->cieY += sd->tLamb.cieY;
746	>	/* gather non-diffuse components */
747	>	i = nr = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
748	>	(rdf != NULL)) ? rdf->ncomp : 0;
749	>	j = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
750	>	(tdf != NULL)) ? tdf->ncomp : 0;
751	>	n = i + j;
752	>	if (n > 0 && (cdarr = (const SDCDst **)malloc(n*sizeof(SDCDst *))) == NULL)
753	>	return SDEmemory;
754	>	while (j-- > 0) {               /* non-diffuse transmission */
755	>	cdarr[i+j] = (*tdf->comp[j].func->getCDist)(inVec, &tdf->comp[j]);
756	>	if (cdarr[i+j] == NULL) {
757	>	free(cdarr);
758	>	return SDEmemory;
759	>	}
760	>	sv->cieY += cdarr[i+j]->cTotal;
761		}
762	<	theta = asin(-xp[2]);
763	<	if (!FEQ(theta,0.0)) {
764	<	*xfp++ = "-ry";
765	<	sprintf(bufs[bn], " %f", theta*(180./PI));
766	<	*xfp++ = bufs[bn++];
762	>	while (i-- > 0) {               /* non-diffuse reflection */
763	>	cdarr[i] = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]);
764	>	if (cdarr[i] == NULL) {
765	>	free(cdarr);
766	>	return SDEmemory;
767	>	}
768	>	sv->cieY += cdarr[i]->cTotal;
769		}
770	<	theta = atan2(xp[1], xp[0]);
771	<	if (!FEQ(theta,0.0)) {
772	<	*xfp++ = "-rz";
773	<	sprintf(bufs[bn], "%f", theta*(180./PI));
660	<	*xfp++ = bufs[bn++];
770	>	if (sv->cieY <= 1e-6) {         /* anything to sample? */
771	>	sv->cieY = .0;
772	>	memset(ioVec, 0, 3*sizeof(double));
773	>	return SDEnone;
774		}
775	<	*xfp = NULL;
776	<	return(xfp - xfarg);
775	>	/* scale random variable */
776	>	randX *= sv->cieY;
777	>	/* diffuse reflection? */
778	>	if (randX < rdiff) {
779	>	SDdiffuseSamp(ioVec, inFront, randX/rdiff);
780	>	goto done;
781	>	}
782	>	randX -= rdiff;
783	>	/* diffuse transmission? */
784	>	if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) {
785	>	if (randX < sd->tLamb.cieY) {
786	>	sv->spec = sd->tLamb.spec;
787	>	SDdiffuseSamp(ioVec, !inFront, randX/sd->tLamb.cieY);
788	>	goto done;
789	>	}
790	>	randX -= sd->tLamb.cieY;
791	>	}
792	>	/* else one of cumulative dist. */
793	>	for (i = 0; i < n && randX < cdarr[i]->cTotal; i++)
794	>	randX -= cdarr[i]->cTotal;
795	>	if (i >= n)
796	>	return SDEinternal;
797	>	/* compute sample direction */
798	>	sdc = (i < nr) ? &rdf->comp[i] : &tdf->comp[i-nr];
799	>	ec = (*sdc->func->sampCDist)(ioVec, randX/cdarr[i]->cTotal, cdarr[i]);
800	>	if (ec)
801	>	return ec;
802	>	/* compute color */
803	>	j = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
804	>	if (j <= 0) {
805	>	sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error",
806	>	sd->name);
807	>	return SDEinternal;
808	>	}
809	>	sv->spec = sdc->cspec[0];
810	>	rdiff = coef[0];
811	>	while (--j) {
812	>	c_cmix(&sv->spec, rdiff, &sv->spec, coef[j], &sdc->cspec[j]);
813	>	rdiff += coef[j];
814	>	}
815	>	done:
816	>	if (cdarr != NULL)
817	>	free(cdarr);
818	>	/* make sure everything is set */
819	>	c_ccvt(&sv->spec, C_CSXY+C_CSSPEC);
820	>	return SDEnone;
821		}
822
823	+	/* Compute World->BSDF transform from surface normal and up (Y) vector */
824	+	SDError
825	+	SDcompXform(RREAL vMtx[3][3], const FVECT sNrm, const FVECT uVec)
826	+	{
827	+	if ((vMtx == NULL) | (sNrm == NULL) | (uVec == NULL))
828	+	return SDEargument;
829	+	VCOPY(vMtx[2], sNrm);
830	+	if (normalize(vMtx[2]) == 0)
831	+	return SDEargument;
832	+	fcross(vMtx[0], uVec, vMtx[2]);
833	+	if (normalize(vMtx[0]) == 0)
834	+	return SDEargument;
835	+	fcross(vMtx[1], vMtx[2], vMtx[0]);
836	+	return SDEnone;
837	+	}
838
839	<	int
840	<	getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
841	<	MAT4 xm,
670	<	FVECT nrm,
671	<	UpDir ud,
672	<	char *xfbuf
673	<	)
839	>	/* Compute inverse transform */
840	>	SDError
841	>	SDinvXform(RREAL iMtx[3][3], RREAL vMtx[3][3])
842		{
843	<	char    *xfargs[7];
844	<	XF      myxf;
677	<	FVECT   updir, xdest, ydest;
678	<	int     i;
843	>	RREAL   mTmp[3][3];
844	>	double  d;
845
846	<	updir[0] = updir[1] = updir[2] = 0.;
847	<	switch (ud) {
848	<	case UDzneg:
849	<	updir[2] = -1.;
850	<	break;
851	<	case UDyneg:
852	<	updir[1] = -1.;
853	<	break;
854	<	case UDxneg:
855	<	updir[0] = -1.;
690	<	break;
691	<	case UDxpos:
692	<	updir[0] = 1.;
693	<	break;
694	<	case UDypos:
695	<	updir[1] = 1.;
696	<	break;
697	<	case UDzpos:
698	<	updir[2] = 1.;
699	<	break;
700	<	case UDunknown:
701	<	return(0);
846	>	if ((iMtx == NULL) | (vMtx == NULL))
847	>	return SDEargument;
848	>	/* compute determinant */
849	>	mTmp[0][0] = vMtx[2][2]*vMtx[1][1] - vMtx[2][1]*vMtx[1][2];
850	>	mTmp[0][1] = vMtx[2][1]*vMtx[0][2] - vMtx[2][2]*vMtx[0][1];
851	>	mTmp[0][2] = vMtx[1][2]*vMtx[0][1] - vMtx[1][1]*vMtx[0][2];
852	>	d = vMtx[0][0]*mTmp[0][0] + vMtx[1][0]*mTmp[0][1] + vMtx[2][0]*mTmp[0][2];
853	>	if (d == 0) {
854	>	strcpy(SDerrorDetail, "Zero determinant in matrix inversion");
855	>	return SDEargument;
856		}
857	<	fcross(xdest, updir, nrm);
858	<	if (normalize(xdest) == 0.0)
859	<	return(0);
860	<	fcross(ydest, nrm, xdest);
861	<	xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
862	<	copymat4(xm, myxf.xfm);
863	<	if (xfbuf == NULL)
864	<	return(1);
865	<	/* return xf arguments as well */
866	<	for (i = 0; xfargs[i] != NULL; i++) {
867	<	*xfbuf++ = ' ';
868	<	strcpy(xfbuf, xfargs[i]);
869	<	while (*xfbuf) ++xfbuf;
857	>	d = 1./d;                       /* invert matrix */
858	>	mTmp[0][0] *= d; mTmp[0][1] *= d; mTmp[0][2] *= d;
859	>	mTmp[1][0] = d*(vMtx[2][0]*vMtx[1][2] - vMtx[2][2]*vMtx[1][0]);
860	>	mTmp[1][1] = d*(vMtx[2][2]*vMtx[0][0] - vMtx[2][0]*vMtx[0][2]);
861	>	mTmp[1][2] = d*(vMtx[1][0]*vMtx[0][2] - vMtx[1][2]*vMtx[0][0]);
862	>	mTmp[2][0] = d*(vMtx[2][1]*vMtx[1][0] - vMtx[2][0]*vMtx[1][1]);
863	>	mTmp[2][1] = d*(vMtx[2][0]*vMtx[0][1] - vMtx[2][1]*vMtx[0][0]);
864	>	mTmp[2][2] = d*(vMtx[1][1]*vMtx[0][0] - vMtx[1][0]*vMtx[0][1]);
865	>	memcpy(iMtx, mTmp, sizeof(mTmp));
866	>	return SDEnone;
867	>	}
868	>
869	>	/* Transform and normalize direction (column) vector */
870	>	SDError
871	>	SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT inpVec)
872	>	{
873	>	FVECT   vTmp;
874	>
875	>	if ((resVec == NULL) | (inpVec == NULL))
876	>	return SDEargument;
877	>	if (vMtx == NULL) {             /* assume they just want to normalize */
878	>	if (resVec != inpVec)
879	>	VCOPY(resVec, inpVec);
880	>	return (normalize(resVec) > 0) ? SDEnone : SDEargument;
881		}
882	<	return(1);
882	>	vTmp[0] = DOT(vMtx[0], inpVec);
883	>	vTmp[1] = DOT(vMtx[1], inpVec);
884	>	vTmp[2] = DOT(vMtx[2], inpVec);
885	>	if (normalize(vTmp) == 0)
886	>	return SDEargument;
887	>	VCOPY(resVec, vTmp);
888	>	return SDEnone;
889		}

Diff Legend

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