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.8 by greg, Sun Sep 26 15:43:26 2010 UTC vs.
Revision 2.48 by greg, Mon Mar 24 04:00:45 2014 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-7 && (b)-(a) <= 1e-7)
50	>	/* Retain BSDFs in cache list */
51	>	int                     SDretainSet = SDretainNone;
52
53	<	// returns the name of the given tag
54	<	#ifdef ezxml_name
55	<	#undef ezxml_name
69	<	static char *
70	<	ezxml_name(ezxml_t xml)
53	>	/* Report any error to the indicated stream */
54	>	SDError
55	>	SDreportError(SDError ec, FILE *fp)
56		{
57	<	if (xml == NULL)
58	<	return(NULL);
59	<	return(xml->name);
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		}
76	–	#endif
75
76	<	// returns the given tag's character content or empty string if none
77	<	#ifdef ezxml_txt
78	<	#undef ezxml_txt
79	<	static char *
82	<	ezxml_txt(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("");
83	<	return(xml->txt);
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		}
88	–	#endif
90
91	<
92	<	static int
93	<	ab_getvec(              /* get vector for this angle basis index */
93	<	FVECT v,
94	<	int ndx,
95	<	void *p
96	<	)
91	>	/* Load geometric dimensions and description (if any) */
92	>	static SDError
93	>	SDloadGeometry(SDData *sd, ezxml_t wtl)
94		{
95	<	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
96	<	int     li;
97	<	double  pol, azi, d;
98	<
99	<	if ((ndx < 0) | (ndx >= ab->nangles))
100	<	return(0);
101	<	for (li = 0; ndx >= ab->lat[li].nphis; li++)
102	<	ndx -= ab->lat[li].nphis;
103	<	pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
104	<	azi = 2.*PI*ndx/ab->lat[li].nphis;
105	<	v[2] = d = cos(pol);
106	<	d = sqrt(1. - d*d);     /* sin(pol) */
107	<	v[0] = cos(azi)*d;
108	<	v[1] = sin(azi)*d;
109	<	return(1);
95	>	ezxml_t         node, matl, geom;
96	>	double          cfact;
97	>	const char      *fmt = NULL, *mgfstr;
98	>
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 */
118	<	FVECT v,
119	<	void *p
120	<	)
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;
124	<	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 */
145	<	int ndx,
146	<	void *p
147	<	)
243	>	/* Allocate new spectral distribution function */
244	>	SDSpectralDF *
245	>	SDnewSpectralDF(int nc)
246		{
247	<	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
150	<	int     li;
151	<	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++)
156	<	ndx -= ab->lat[li].nphis;
157	<	theta1 = PI/180. * ab->lat[li+1].tmin;
158	<	if (ab->lat[li].nphis == 1) {           /* special case */
159	<	if (ab->lat[li].tmin > FTINY)
160	<	error(USER, "unsupported BSDF coordinate system");
161	<	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 */
170	<	FVECT v,
171	<	int ndx,
172	<	void *p
173	<	)
267	>	/* Add component(s) to spectral distribution function */
268	>	SDSpectralDF *
269	>	SDaddComponent(SDSpectralDF *odf, int nadd)
270		{
271	<	if (!ab_getvec(v, ndx, p))
176	<	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 */
188	<	FVECT v,
189	<	void *p
190	<	)
290	>	/* Free cached cumulative distributions for BSDF component */
291	>	void
292	>	SDfreeCumulativeCache(SDSpectralDF *df)
293		{
294	<	FVECT  v2;
295	<
194	<	v2[0] = -v[0];
195	<	v2[1] = -v[1];
196	<	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 */
204	<	ezxml_t wab
205	<	)
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;
210	<
211	<	if (!abname || !*abname)
310	>	int     n;
311	>
312	>	if (df == NULL)
313		return;
314	<	for (i = nabases; i--; )
315	<	if (!strcmp(abname, abase_list[i].name))
316	<	return;         /* assume it's the same */
317	<	if (nabases >= MAXABASES)
318	<	error(INTERNAL, "too many angle bases");
218	<	strcpy(abase_list[nabases].name, abname);
219	<	abase_list[nabases].nangles = 0;
220	<	for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
221	<	wbb != NULL; i++, wbb = wbb->next) {
222	<	if (i >= MAXLATS)
223	<	error(INTERNAL, "too many latitudes in custom basis");
224	<	abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
225	<	ezxml_child(ezxml_child(wbb,
226	<	"ThetaBounds"), "UpperTheta")));
227	<	if (!i)
228	<	abase_list[nabases].lat[i].tmin =
229	<	-abase_list[nabases].lat[i+1].tmin;
230	<	else if (!FEQ(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
231	<	"ThetaBounds"), "LowerTheta"))),
232	<	abase_list[nabases].lat[i].tmin))
233	<	error(WARNING, "theta values disagree in custom basis");
234	<	abase_list[nabases].nangles +=
235	<	abase_list[nabases].lat[i].nphis =
236	<	atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
237	<	}
238	<	abase_list[nabases++].lat[i].nphis = 0;
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 */
244	<	const char *unit
245	<	)
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) */
260	<	struct BSDF_data *dp,
261	<	ezxml_t wdb
262	<	)
339	>	/* Initialize an unused BSDF struct (simply clears to zeroes) */
340	>	void
341	>	SDclearBSDF(SDData *sd, const char *fname)
342		{
343	<	ezxml_t         geom;
265	<	double          cfact;
266	<	const char      *fmt, *mgfstr;
267	<
268	<	dp->dim[0] = dp->dim[1] = dp->dim[2] = 0;
269	<	dp->mgf = NULL;
270	<	if ((geom = ezxml_child(wdb, "Width")) != NULL)
271	<	dp->dim[0] = atof(ezxml_txt(geom)) *
272	<	to_meters(ezxml_attr(geom, "unit"));
273	<	if ((geom = ezxml_child(wdb, "Height")) != NULL)
274	<	dp->dim[1] = atof(ezxml_txt(geom)) *
275	<	to_meters(ezxml_attr(geom, "unit"));
276	<	if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
277	<	dp->dim[2] = atof(ezxml_txt(geom)) *
278	<	to_meters(ezxml_attr(geom, "unit"));
279	<	if ((geom = ezxml_child(wdb, "Geometry")) == NULL ||
280	<	(mgfstr = ezxml_txt(geom)) == NULL)
343	>	if (sd == NULL)
344		return;
345	<	if ((fmt = ezxml_attr(geom, "format")) != NULL &&
346	<	strcasecmp(fmt, "MGF")) {
284	<	sprintf(errmsg, "unrecognized geometry format '%s'", fmt);
285	<	error(WARNING, errmsg);
345	>	memset(sd, 0, sizeof(SDData));
346	>	if (fname == NULL)
347		return;
348	<	}
288	<	cfact = to_meters(ezxml_attr(geom, "unit"));
289	<	dp->mgf = (char *)malloc(strlen(mgfstr)+32);
290	<	if (dp->mgf == NULL)
291	<	error(SYSTEM, "out of memory in load_geometry");
292	<	if (cfact < 0.99 || cfact > 1.01)
293	<	sprintf(dp->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
294	<	else
295	<	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 */
301	<	struct BSDF_data *dp,
302	<	ezxml_t wdb
303	<	)
351	>	/* Free data associated with BSDF struct */
352	>	void
353	>	SDfreeBSDF(SDData *sd)
354		{
355	<	char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
306	<	char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
307	<	char  *sdata;
308	<	int  i;
309	<
310	<	if ((!cbasis || !*cbasis) | (!rbasis || !*rbasis)) {
311	<	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 (!strcmp(cbasis, abase_list[i].name)) {
363	<	dp->ninc = abase_list[i].nangles;
317	<	dp->ib_priv = (void *)&abase_list[i];
318	<	dp->ib_vec = ab_getvecR;
319	<	dp->ib_ndx = ab_getndxR;
320	<	dp->ib_ohm = ab_getohm;
321	<	break;
322	<	}
323	<	if (i < 0) {
324	<	sprintf(errmsg, "undefined ColumnAngleBasis '%s'", cbasis);
325	<	error(WARNING, errmsg);
326	<	return;
361	>	if (sd->rf != NULL) {
362	>	SDfreeSpectralDF(sd->rf);
363	>	sd->rf = NULL;
364		}
365	<	for (i = nabases; i--; )
366	<	if (!strcmp(rbasis, abase_list[i].name)) {
367	<	dp->nout = abase_list[i].nangles;
331	<	dp->ob_priv = (void *)&abase_list[i];
332	<	dp->ob_vec = ab_getvec;
333	<	dp->ob_ndx = ab_getndx;
334	<	dp->ob_ohm = ab_getohm;
335	<	break;
336	<	}
337	<	if (i < 0) {
338	<	sprintf(errmsg, "undefined RowAngleBasis '%s'", cbasis);
339	<	error(WARNING, errmsg);
340	<	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) {
345	<	error(WARNING, "missing BSDF ScatteringData");
346	<	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");
351	<	for (i = 0; i < dp->ninc*dp->nout; i++) {
352	<	char  *sdnext = fskip(sdata);
353	<	if (sdnext == NULL) {
354	<	error(WARNING, "bad/missing BSDF ScatteringData");
355	<	free(dp->bsdf); dp->bsdf = NULL;
356	<	return;
357	<	}
358	<	while (*sdnext && isspace(*sdnext))
359	<	sdnext++;
360	<	if (*sdnext == ',') sdnext++;
361	<	dp->bsdf[i] = atof(sdata);
362	<	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);
370	<	}
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 */
376	<	struct BSDF_data *dp
377	<	)
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;
381	<	int             nneg;
382	<	FVECT           v;
383	<	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;
393	<	for (i = dp->ninc; i--; ) {
394	<	dom = getBSDF_incohm(dp,i);
395	<	if (dom <= .0) {
396	<	error(WARNING, "zero/negative incoming solid angle");
397	<	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	>	if ((sd = SDgetCache(fname)) == NULL) {
426	>	SDreportError(SDEmemory, stderr);
427	>	return NULL;
428		}
429	<	if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
430	<	sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
431	<	100.*(hemi_total/PI - 1.));
432	<	error(WARNING, errmsg);
429	>	if (!SDisLoaded(sd) && (ec = SDloadFile(sd, fname))) {
430	>	SDreportError(ec, stderr);
431	>	SDfreeCache(sd);
432	>	return NULL;
433		}
434	<	dom = PI / hemi_total;          /* fix normalization */
435	<	for (i = dp->ninc; i--; )
413	<	omega_iarr[i] *= dom;
414	<	/* outgoing projected solid angles */
415	<	hemi_total = .0;
416	<	for (o = dp->nout; o--; ) {
417	<	dom = getBSDF_outohm(dp,o);
418	<	if (dom <= .0) {
419	<	error(WARNING, "zero/negative outgoing solid angle");
420	<	continue;
421	<	}
422	<	if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
423	<	error(WARNING, "illegal outgoing BSDF direction");
424	<	free(omega_iarr); free(omega_oarr);
425	<	return(0);
426	<	}
427	<	hemi_total += omega_oarr[o] = dom * v[2];
428	<	}
429	<	if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
430	<	sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
431	<	100.*(hemi_total/PI - 1.));
432	<	error(WARNING, errmsg);
433	<	}
434	<	dom = PI / hemi_total;          /* fix normalization */
435	<	for (o = dp->nout; o--; )
436	<	omega_oarr[o] *= dom;
437	<	nneg = 0;                       /* check outgoing totals */
438	<	for (i = 0; i < dp->ninc; i++) {
439	<	hemi_total = .0;
440	<	for (o = dp->nout; o--; ) {
441	<	double  f = BSDF_value(dp,i,o);
442	<	if (f >= .0)
443	<	hemi_total += f*omega_oarr[o];
444	<	else {
445	<	nneg += (f < -FTINY);
446	<	BSDF_value(dp,i,o) = .0f;
447	<	}
448	<	}
449	<	if (hemi_total > 1.01) {
450	<	sprintf(errmsg,
451	<	"incoming BSDF direction %d passes %.1f%% of light",
452	<	i, 100.*hemi_total);
453	<	error(WARNING, errmsg);
454	<	}
455	<	}
456	<	if (nneg) {
457	<	sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
458	<	error(WARNING, errmsg);
459	<	}
460	<	full_total = .0;                /* reverse roles and check again */
461	<	for (o = 0; o < dp->nout; o++) {
462	<	hemi_total = .0;
463	<	for (i = dp->ninc; i--; )
464	<	hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
434	>	return sd;
435	>	}
436
437	<	if (hemi_total > 1.01) {
438	<	sprintf(errmsg,
439	<	"outgoing BSDF direction %d collects %.1f%% of light",
440	<	o, 100.*hemi_total);
441	<	error(WARNING, errmsg);
437	>	/* Free a BSDF from our cache (clear all if NULL) */
438	>	void
439	>	SDfreeCache(const SDData *sd)
440	>	{
441	>	struct SDCache_s        *sdl, *sdLast = NULL;
442	>
443	>	if (sd == NULL) {               /* free entire list */
444	>	while ((sdl = SDcacheList) != NULL) {
445	>	SDcacheList = sdl->next;
446	>	SDfreeBSDF(&sdl->bsdf);
447	>	free(sdl);
448		}
449	<	full_total += hemi_total*omega_oarr[o];
449	>	return;
450		}
451	<	full_total /= PI;
452	<	if (full_total > 1.00001) {
453	<	sprintf(errmsg, "BSDF transfers %.4f%% of light",
454	<	100.*full_total);
455	<	error(WARNING, errmsg);
451	>	for (sdl = SDcacheList; sdl != NULL; sdl = (sdLast=sdl)->next)
452	>	if (&sdl->bsdf == sd)
453	>	break;
454	>	if (sdl == NULL || (sdl->refcnt -= (sdl->refcnt > 0)))
455	>	return;                 /* missing or still in use */
456	>	/* keep unreferenced data? */
457	>	if (SDisLoaded(sd) && SDretainSet) {
458	>	if (SDretainSet == SDretainAll)
459	>	return;         /* keep everything */
460	>	/* else free cumulative data */
461	>	SDfreeCumulativeCache(sd->rf);
462	>	SDfreeCumulativeCache(sd->rb);
463	>	SDfreeCumulativeCache(sd->tf);
464	>	SDfreeCumulativeCache(sd->tb);
465	>	return;
466		}
467	<	free(omega_iarr); free(omega_oarr);
468	<	return(1);
467	>	/* remove from list and free */
468	>	if (sdLast == NULL)
469	>	SDcacheList = sdl->next;
470	>	else
471	>	sdLast->next = sdl->next;
472	>	SDfreeBSDF(&sdl->bsdf);
473	>	free(sdl);
474		}
475
476	<
477	<	struct BSDF_data *
478	<	load_BSDF(              /* load BSDF data from file */
487	<	char *fname
488	<	)
476	>	/* Sample an individual BSDF component */
477	>	SDError
478	>	SDsampComponent(SDValue *sv, FVECT ioVec, double randX, SDComponent *sdc)
479		{
480	<	char                    *path;
481	<	ezxml_t                 fl, wtl, wld, wdb;
482	<	struct BSDF_data        *dp;
483	<
484	<	path = getpath(fname, getrlibpath(), R_OK);
485	<	if (path == NULL) {
486	<	sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
487	<	error(WARNING, errmsg);
488	<	return(NULL);
480	>	float           coef[SDmaxCh];
481	>	SDError         ec;
482	>	FVECT           inVec;
483	>	const SDCDst    *cd;
484	>	double          d;
485	>	int             n;
486	>	/* check arguments */
487	>	if ((sv == NULL) | (ioVec == NULL) | (sdc == NULL))
488	>	return SDEargument;
489	>	/* get cumulative distribution */
490	>	VCOPY(inVec, ioVec);
491	>	sv->cieY = 0;
492	>	cd = (*sdc->func->getCDist)(inVec, sdc);
493	>	if (cd != NULL)
494	>	sv->cieY = cd->cTotal;
495	>	if (sv->cieY <= 1e-6) {         /* nothing to sample? */
496	>	sv->spec = c_dfcolor;
497	>	memset(ioVec, 0, 3*sizeof(double));
498	>	return SDEnone;
499		}
500	<	fl = ezxml_parse_file(path);
501	<	if (fl == NULL) {
502	<	sprintf(errmsg, "cannot open BSDF \"%s\"", path);
503	<	error(WARNING, errmsg);
504	<	return(NULL);
500	>	/* compute sample direction */
501	>	ec = (*sdc->func->sampCDist)(ioVec, randX, cd);
502	>	if (ec)
503	>	return ec;
504	>	/* get BSDF color */
505	>	n = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
506	>	if (n <= 0) {
507	>	strcpy(SDerrorDetail, "BSDF sample value error");
508	>	return SDEinternal;
509		}
510	<	if (ezxml_error(fl)[0]) {
511	<	sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
512	<	error(WARNING, errmsg);
513	<	ezxml_free(fl);
514	<	return(NULL);
510	>	sv->spec = sdc->cspec[0];
511	>	d = coef[0];
512	>	while (--n) {
513	>	c_cmix(&sv->spec, d, &sv->spec, coef[n], &sdc->cspec[n]);
514	>	d += coef[n];
515		}
516	<	if (strcmp(ezxml_name(fl), "WindowElement")) {
517	<	sprintf(errmsg,
518	<	"BSDF \"%s\": top level node not 'WindowElement'",
515	<	path);
516	<	error(WARNING, errmsg);
517	<	ezxml_free(fl);
518	<	return(NULL);
519	<	}
520	<	wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
521	<	load_angle_basis(ezxml_child(ezxml_child(wtl,
522	<	"DataDefinition"), "AngleBasis"));
523	<	dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
524	<	load_geometry(dp, ezxml_child(wtl, "Material"));
525	<	for (wld = ezxml_child(wtl, "WavelengthData");
526	<	wld != NULL; wld = wld->next) {
527	<	if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
528	<	continue;
529	<	wdb = ezxml_child(wld, "WavelengthDataBlock");
530	<	if (wdb == NULL) continue;
531	<	if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
532	<	"Transmission Front"))
533	<	continue;
534	<	load_bsdf_data(dp, wdb);        /* load front BTDF */
535	<	break;                          /* ignore the rest */
536	<	}
537	<	ezxml_free(fl);                         /* done with XML file */
538	<	if (!check_bsdf_data(dp)) {
539	<	sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
540	<	error(WARNING, errmsg);
541	<	free_BSDF(dp);
542	<	dp = NULL;
543	<	}
544	<	return(dp);
516	>	/* make sure everything is set */
517	>	c_ccvt(&sv->spec, C_CSXY+C_CSSPEC);
518	>	return SDEnone;
519		}
520
521	+	#define MS_MAXDIM       15
522
523	+	/* Convert 1-dimensional random variable to N-dimensional */
524		void
525	<	free_BSDF(              /* free BSDF data structure */
550	<	struct BSDF_data *b
551	<	)
525	>	SDmultiSamp(double t[], int n, double randX)
526		{
527	<	if (b == NULL)
527	>	unsigned        nBits;
528	>	double          scale;
529	>	bitmask_t       ndx, coord[MS_MAXDIM];
530	>
531	>	if (n <= 0)                     /* check corner cases */
532		return;
533	<	if (b->mgf != NULL)
534	<	free(b->mgf);
535	<	if (b->bsdf != NULL)
536	<	free(b->bsdf);
537	<	free(b);
533	>	if (randX < 0) randX = 0;
534	>	else if (randX >= 1.) randX = 0.999999999999999;
535	>	if (n == 1) {
536	>	t[0] = randX;
537	>	return;
538	>	}
539	>	while (n > MS_MAXDIM)           /* punt for higher dimensions */
540	>	t[--n] = rand()*(1./(RAND_MAX+.5));
541	>	nBits = (8*sizeof(bitmask_t) - 1) / n;
542	>	ndx = randX * (double)((bitmask_t)1 << (nBits*n));
543	>	/* get coordinate on Hilbert curve */
544	>	hilbert_i2c(n, nBits, ndx, coord);
545	>	/* convert back to [0,1) range */
546	>	scale = 1. / (double)((bitmask_t)1 << nBits);
547	>	while (n--)
548	>	t[n] = scale * ((double)coord[n] + rand()*(1./(RAND_MAX+.5)));
549		}
550
551	+	#undef MS_MAXDIM
552
553	<	int
554	<	r_BSDF_incvec(          /* compute random input vector at given location */
555	<	FVECT v,
566	<	struct BSDF_data *b,
567	<	int i,
568	<	double rv,
569	<	MAT4 xm
570	<	)
553	>	/* Generate diffuse hemispherical sample */
554	>	static void
555	>	SDdiffuseSamp(FVECT outVec, int outFront, double randX)
556		{
557	<	FVECT   pert;
558	<	double  rad;
559	<	int     j;
560	<
561	<	if (!getBSDF_incvec(v, b, i))
562	<	return(0);
563	<	rad = sqrt(getBSDF_incohm(b, i) / PI);
579	<	multisamp(pert, 3, rv);
580	<	for (j = 0; j < 3; j++)
581	<	v[j] += rad*(2.*pert[j] - 1.);
582	<	if (xm != NULL)
583	<	multv3(v, v, xm);
584	<	return(normalize(v) != 0.0);
557	>	/* convert to position on hemisphere */
558	>	SDmultiSamp(outVec, 2, randX);
559	>	SDsquare2disk(outVec, outVec[0], outVec[1]);
560	>	outVec[2] = 1. - outVec[0]*outVec[0] - outVec[1]*outVec[1];
561	>	outVec[2] = sqrt(outVec[2]*(outVec[2]>0));
562	>	if (!outFront)                  /* going out back? */
563	>	outVec[2] = -outVec[2];
564		}
565
566	+	/* Query projected solid angle coverage for non-diffuse BSDF direction */
567	+	SDError
568	+	SDsizeBSDF(double *projSA, const FVECT v1, const RREAL *v2,
569	+	int qflags, const SDData *sd)
570	+	{
571	+	SDSpectralDF    *rdf, *tdf;
572	+	SDError         ec;
573	+	int             i;
574	+	/* check arguments */
575	+	if ((projSA == NULL) | (v1 == NULL) | (sd == NULL))
576	+	return SDEargument;
577	+	/* initialize extrema */
578	+	switch (qflags) {
579	+	case SDqueryMax:
580	+	projSA[0] = .0;
581	+	break;
582	+	case SDqueryMin+SDqueryMax:
583	+	projSA[1] = .0;
584	+	/* fall through */
585	+	case SDqueryMin:
586	+	projSA[0] = 10.;
587	+	break;
588	+	case 0:
589	+	return SDEargument;
590	+	}
591	+	if (v1[2] > 0) {                /* front surface query? */
592	+	rdf = sd->rf;
593	+	tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
594	+	} else {
595	+	rdf = sd->rb;
596	+	tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
597	+	}
598	+	if (v2 != NULL)                 /* bidirectional? */
599	+	if (v1[2] > 0 ^ v2[2] > 0)
600	+	rdf = NULL;
601	+	else
602	+	tdf = NULL;
603	+	ec = SDEdata;                   /* run through components */
604	+	for (i = (rdf==NULL) ? 0 : rdf->ncomp; i--; ) {
605	+	ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2,
606	+	qflags, &rdf->comp[i]);
607	+	if (ec)
608	+	return ec;
609	+	}
610	+	for (i = (tdf==NULL) ? 0 : tdf->ncomp; i--; ) {
611	+	ec = (*tdf->comp[i].func->queryProjSA)(projSA, v1, v2,
612	+	qflags, &tdf->comp[i]);
613	+	if (ec)
614	+	return ec;
615	+	}
616	+	if (ec) {                       /* all diffuse? */
617	+	projSA[0] = M_PI;
618	+	if (qflags == SDqueryMin+SDqueryMax)
619	+	projSA[1] = M_PI;
620	+	} else if (qflags == SDqueryMin+SDqueryMax && projSA[0] > projSA[1])
621	+	projSA[0] = projSA[1];
622	+	return SDEnone;
623	+	}
624
625	<	int
626	<	r_BSDF_outvec(          /* compute random output vector at given location */
627	<	FVECT v,
591	<	struct BSDF_data *b,
592	<	int o,
593	<	double rv,
594	<	MAT4 xm
595	<	)
625	>	/* Return BSDF for the given incident and scattered ray vectors */
626	>	SDError
627	>	SDevalBSDF(SDValue *sv, const FVECT outVec, const FVECT inVec, const SDData *sd)
628		{
629	<	FVECT   pert;
630	<	double  rad;
631	<	int     j;
632	<
633	<	if (!getBSDF_outvec(v, b, o))
634	<	return(0);
635	<	rad = sqrt(getBSDF_outohm(b, o) / PI);
636	<	multisamp(pert, 3, rv);
637	<	for (j = 0; j < 3; j++)
638	<	v[j] += rad*(2.*pert[j] - 1.);
639	<	if (xm != NULL)
640	<	multv3(v, v, xm);
641	<	return(normalize(v) != 0.0);
629	>	int             inFront, outFront;
630	>	SDSpectralDF    *sdf;
631	>	float           coef[SDmaxCh];
632	>	int             nch, i;
633	>	/* check arguments */
634	>	if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL))
635	>	return SDEargument;
636	>	/* whose side are we on? */
637	>	inFront = (inVec[2] > 0);
638	>	outFront = (outVec[2] > 0);
639	>	/* start with diffuse portion */
640	>	if (inFront & outFront) {
641	>	*sv = sd->rLambFront;
642	>	sdf = sd->rf;
643	>	} else if (!(inFront | outFront)) {
644	>	*sv = sd->rLambBack;
645	>	sdf = sd->rb;
646	>	} else if (inFront) {
647	>	*sv = sd->tLamb;
648	>	sdf = (sd->tf != NULL) ? sd->tf : sd->tb;
649	>	} else /* inBack */ {
650	>	*sv = sd->tLamb;
651	>	sdf = (sd->tb != NULL) ? sd->tb : sd->tf;
652	>	}
653	>	sv->cieY *= 1./M_PI;
654	>	/* add non-diffuse components */
655	>	i = (sdf != NULL) ? sdf->ncomp : 0;
656	>	while (i-- > 0) {
657	>	nch = (*sdf->comp[i].func->getBSDFs)(coef, outVec, inVec,
658	>	&sdf->comp[i]);
659	>	while (nch-- > 0) {
660	>	c_cmix(&sv->spec, sv->cieY, &sv->spec,
661	>	coef[nch], &sdf->comp[i].cspec[nch]);
662	>	sv->cieY += coef[nch];
663	>	}
664	>	}
665	>	/* make sure everything is set */
666	>	c_ccvt(&sv->spec, C_CSXY+C_CSSPEC);
667	>	return SDEnone;
668		}
669
670	<
671	<	static int
672	<	addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
615	<	char *xfarg[],
616	<	FVECT xp,
617	<	FVECT yp,
618	<	FVECT zp
619	<	)
670	>	/* Compute directional hemispherical scattering at this incident angle */
671	>	double
672	>	SDdirectHemi(const FVECT inVec, int sflags, const SDData *sd)
673		{
674	<	static char     bufs[3][16];
675	<	int     bn = 0;
676	<	char    **xfp = xfarg;
677	<	double  theta;
674	>	double          hsum;
675	>	SDSpectralDF    *rdf, *tdf;
676	>	const SDCDst    *cd;
677	>	int             i;
678	>	/* check arguments */
679	>	if ((inVec == NULL) | (sd == NULL))
680	>	return .0;
681	>	/* gather diffuse components */
682	>	if (inVec[2] > 0) {
683	>	hsum = sd->rLambFront.cieY;
684	>	rdf = sd->rf;
685	>	tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
686	>	} else /* !inFront */ {
687	>	hsum = sd->rLambBack.cieY;
688	>	rdf = sd->rb;
689	>	tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
690	>	}
691	>	if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
692	>	hsum = .0;
693	>	if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
694	>	hsum += sd->tLamb.cieY;
695	>	/* gather non-diffuse components */
696	>	i = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
697	>	(rdf != NULL)) ? rdf->ncomp : 0;
698	>	while (i-- > 0) {               /* non-diffuse reflection */
699	>	cd = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]);
700	>	if (cd != NULL)
701	>	hsum += cd->cTotal;
702	>	}
703	>	i = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
704	>	(tdf != NULL)) ? tdf->ncomp : 0;
705	>	while (i-- > 0) {               /* non-diffuse transmission */
706	>	cd = (*tdf->comp[i].func->getCDist)(inVec, &tdf->comp[i]);
707	>	if (cd != NULL)
708	>	hsum += cd->cTotal;
709	>	}
710	>	return hsum;
711	>	}
712
713	<	if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
714	<	/* Special case for X' along Z-axis */
715	<	theta = -atan2(yp[0], yp[1]);
716	<	*xfp++ = "-ry";
717	<	*xfp++ = xp[2] < 0.0 ? "90" : "-90";
718	<	*xfp++ = "-rz";
719	<	sprintf(bufs[bn], "%f", theta*(180./PI));
720	<	*xfp++ = bufs[bn++];
721	<	return(xfp - xfarg);
713	>	/* Sample BSDF direction based on the given random variable */
714	>	SDError
715	>	SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int sflags, const SDData *sd)
716	>	{
717	>	SDError         ec;
718	>	FVECT           inVec;
719	>	int             inFront;
720	>	SDSpectralDF    *rdf, *tdf;
721	>	double          rdiff;
722	>	float           coef[SDmaxCh];
723	>	int             i, j, n, nr;
724	>	SDComponent     *sdc;
725	>	const SDCDst    **cdarr = NULL;
726	>	/* check arguments */
727	>	if ((sv == NULL) | (ioVec == NULL) | (sd == NULL) |
728	>	(randX < 0) | (randX >= 1.))
729	>	return SDEargument;
730	>	/* whose side are we on? */
731	>	VCOPY(inVec, ioVec);
732	>	inFront = (inVec[2] > 0);
733	>	/* remember diffuse portions */
734	>	if (inFront) {
735	>	*sv = sd->rLambFront;
736	>	rdf = sd->rf;
737	>	tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
738	>	} else /* !inFront */ {
739	>	*sv = sd->rLambBack;
740	>	rdf = sd->rb;
741	>	tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
742		}
743	<	theta = atan2(yp[2], zp[2]);
744	<	if (!FEQ(theta,0.0)) {
745	<	*xfp++ = "-rx";
746	<	sprintf(bufs[bn], "%f", theta*(180./PI));
747	<	*xfp++ = bufs[bn++];
743	>	if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
744	>	sv->cieY = .0;
745	>	rdiff = sv->cieY;
746	>	if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
747	>	sv->cieY += sd->tLamb.cieY;
748	>	/* gather non-diffuse components */
749	>	i = nr = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
750	>	(rdf != NULL)) ? rdf->ncomp : 0;
751	>	j = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
752	>	(tdf != NULL)) ? tdf->ncomp : 0;
753	>	n = i + j;
754	>	if (n > 0 && (cdarr = (const SDCDst **)malloc(n*sizeof(SDCDst *))) == NULL)
755	>	return SDEmemory;
756	>	while (j-- > 0) {               /* non-diffuse transmission */
757	>	cdarr[i+j] = (*tdf->comp[j].func->getCDist)(inVec, &tdf->comp[j]);
758	>	if (cdarr[i+j] == NULL)
759	>	cdarr[i+j] = &SDemptyCD;
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	>	cdarr[i] = &SDemptyCD;
766	>	else
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));
652	<	*xfp++ = bufs[bn++];
769	>	if (sv->cieY <= 1e-6) {         /* anything to sample? */
770	>	sv->cieY = .0;
771	>	memset(ioVec, 0, 3*sizeof(double));
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	>	/* make sure everything is set */
818	>	c_ccvt(&sv->spec, C_CSXY+C_CSSPEC);
819	>	return SDEnone;
820		}
821
822	+	/* Compute World->BSDF transform from surface normal and up (Y) vector */
823	+	SDError
824	+	SDcompXform(RREAL vMtx[3][3], const FVECT sNrm, const FVECT uVec)
825	+	{
826	+	if ((vMtx == NULL) | (sNrm == NULL) | (uVec == NULL))
827	+	return SDEargument;
828	+	VCOPY(vMtx[2], sNrm);
829	+	if (normalize(vMtx[2]) == 0)
830	+	return SDEargument;
831	+	fcross(vMtx[0], uVec, vMtx[2]);
832	+	if (normalize(vMtx[0]) == 0)
833	+	return SDEargument;
834	+	fcross(vMtx[1], vMtx[2], vMtx[0]);
835	+	return SDEnone;
836	+	}
837
838	<	int
839	<	getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
840	<	MAT4 xm,
662	<	FVECT nrm,
663	<	UpDir ud,
664	<	char *xfbuf
665	<	)
838	>	/* Compute inverse transform */
839	>	SDError
840	>	SDinvXform(RREAL iMtx[3][3], RREAL vMtx[3][3])
841		{
842	<	char    *xfargs[7];
843	<	XF      myxf;
669	<	FVECT   updir, xdest, ydest;
670	<	int     i;
842	>	RREAL   mTmp[3][3];
843	>	double  d;
844
845	<	updir[0] = updir[1] = updir[2] = 0.;
846	<	switch (ud) {
847	<	case UDzneg:
848	<	updir[2] = -1.;
849	<	break;
850	<	case UDyneg:
851	<	updir[1] = -1.;
852	<	break;
853	<	case UDxneg:
854	<	updir[0] = -1.;
682	<	break;
683	<	case UDxpos:
684	<	updir[0] = 1.;
685	<	break;
686	<	case UDypos:
687	<	updir[1] = 1.;
688	<	break;
689	<	case UDzpos:
690	<	updir[2] = 1.;
691	<	break;
692	<	case UDunknown:
693	<	return(0);
845	>	if ((iMtx == NULL) | (vMtx == NULL))
846	>	return SDEargument;
847	>	/* compute determinant */
848	>	mTmp[0][0] = vMtx[2][2]*vMtx[1][1] - vMtx[2][1]*vMtx[1][2];
849	>	mTmp[0][1] = vMtx[2][1]*vMtx[0][2] - vMtx[2][2]*vMtx[0][1];
850	>	mTmp[0][2] = vMtx[1][2]*vMtx[0][1] - vMtx[1][1]*vMtx[0][2];
851	>	d = vMtx[0][0]*mTmp[0][0] + vMtx[1][0]*mTmp[0][1] + vMtx[2][0]*mTmp[0][2];
852	>	if (d == 0) {
853	>	strcpy(SDerrorDetail, "Zero determinant in matrix inversion");
854	>	return SDEargument;
855		}
856	<	fcross(xdest, updir, nrm);
857	<	if (normalize(xdest) == 0.0)
858	<	return(0);
859	<	fcross(ydest, nrm, xdest);
860	<	xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
861	<	copymat4(xm, myxf.xfm);
862	<	if (xfbuf == NULL)
863	<	return(1);
864	<	/* return xf arguments as well */
865	<	for (i = 0; xfargs[i] != NULL; i++) {
866	<	*xfbuf++ = ' ';
867	<	strcpy(xfbuf, xfargs[i]);
868	<	while (*xfbuf) ++xfbuf;
856	>	d = 1./d;                       /* invert matrix */
857	>	mTmp[0][0] *= d; mTmp[0][1] *= d; mTmp[0][2] *= d;
858	>	mTmp[1][0] = d*(vMtx[2][0]*vMtx[1][2] - vMtx[2][2]*vMtx[1][0]);
859	>	mTmp[1][1] = d*(vMtx[2][2]*vMtx[0][0] - vMtx[2][0]*vMtx[0][2]);
860	>	mTmp[1][2] = d*(vMtx[1][0]*vMtx[0][2] - vMtx[1][2]*vMtx[0][0]);
861	>	mTmp[2][0] = d*(vMtx[2][1]*vMtx[1][0] - vMtx[2][0]*vMtx[1][1]);
862	>	mTmp[2][1] = d*(vMtx[2][0]*vMtx[0][1] - vMtx[2][1]*vMtx[0][0]);
863	>	mTmp[2][2] = d*(vMtx[1][1]*vMtx[0][0] - vMtx[1][0]*vMtx[0][1]);
864	>	memcpy(iMtx, mTmp, sizeof(mTmp));
865	>	return SDEnone;
866	>	}
867	>
868	>	/* Transform and normalize direction (column) vector */
869	>	SDError
870	>	SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT inpVec)
871	>	{
872	>	FVECT   vTmp;
873	>
874	>	if ((resVec == NULL) | (inpVec == NULL))
875	>	return SDEargument;
876	>	if (vMtx == NULL) {             /* assume they just want to normalize */
877	>	if (resVec != inpVec)
878	>	VCOPY(resVec, inpVec);
879	>	return (normalize(resVec) > 0) ? SDEnone : SDEargument;
880		}
881	<	return(1);
881	>	vTmp[0] = DOT(vMtx[0], inpVec);
882	>	vTmp[1] = DOT(vMtx[1], inpVec);
883	>	vTmp[2] = DOT(vMtx[2], inpVec);
884	>	if (normalize(vTmp) == 0)
885	>	return SDEargument;
886	>	VCOPY(resVec, vTmp);
887	>	return SDEnone;
888		}

Diff Legend

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