ViewVC Help

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

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

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

Diff Legend

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