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

Diff Legend

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