ViewVC Help

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

root/Development/ray/src/common/bsdf_m.c

Comparing ray/src/common/bsdf_m.c (file contents):
Revision 3.6 by greg, Mon Feb 21 22:50:37 2011 UTC vs.
Revision 3.37 by greg, Fri Jan 5 21:00:24 2018 UTC

#	Line 11 | Line 11 | static const char RCSid[] = "$Id$";
11		*
12		*/
13
14	+	#define _USE_MATH_DEFINES
15		#include "rtio.h"
15	–	#include <stdlib.h>
16		#include <math.h>
17		#include <ctype.h>
18		#include "ezxml.h"
#	Line 28 | Line 28 | static const char RCSid[] = "$Id$";
28		#define RC_INTERR       (-4)
29		#define RC_MEMERR       (-5)
30
31	<	#define MAXLATS         46              /* maximum number of latitudes */
32	<
33	<	/* BSDF angle specification */
34	<	typedef struct {
35	<	char    name[64];               /* basis name */
36	<	int     nangles;                /* total number of directions */
37	<	struct {
38	<	float   tmin;                   /* starting theta */
39	<	int     nphis;                  /* number of phis (0 term) */
40	<	} lat[MAXLATS+1];               /* latitudes */
41	<	} ANGLE_BASIS;
42	<
43	<	#define MAXABASES       7               /* limit on defined bases */
44	<
45	<	static ANGLE_BASIS      abase_list[MAXABASES] = {
31	>	ANGLE_BASIS     abase_list[MAXABASES] = {
32		{
33		"LBNL/Klems Full", 145,
34	<	{ {-5., 1},
34	>	{ {0., 1},
35		{5., 8},
36		{15., 16},
37		{25., 20},
#	Line 57 | Line 43 | static ANGLE_BASIS     abase_list[MAXABASES] = {
43		{90., 0} }
44		}, {
45		"LBNL/Klems Half", 73,
46	<	{ {-6.5, 1},
46	>	{ {0., 1},
47		{6.5, 8},
48		{19.5, 12},
49		{32.5, 16},
#	Line 67 | Line 53 | static ANGLE_BASIS     abase_list[MAXABASES] = {
53		{90., 0} }
54		}, {
55		"LBNL/Klems Quarter", 41,
56	<	{ {-9., 1},
56	>	{ {0., 1},
57		{9., 8},
58		{27., 12},
59		{46., 12},
#	Line 76 | Line 62 | static ANGLE_BASIS     abase_list[MAXABASES] = {
62		}
63		};
64
65	<	static int      nabases = 3;    /* current number of defined bases */
65	>	int             nabases = 3;            /* current number of defined bases */
66
67	+	C_COLOR         mtx_RGB_prim[3];        /* our RGB primaries  */
68	+	float           mtx_RGB_coef[3];        /* corresponding Y coefficients */
69	+
70	+	enum {mtx_Y, mtx_X, mtx_Z};             /* matrix components (mtx_Y==0) */
71	+
72	+	/* check if two real values are near enough to equal */
73		static int
74		fequal(double a, double b)
75		{
76	<	if (b != .0)
76	>	if (b != 0)
77		a = a/b - 1.;
78		return (a <= 1e-6) & (a >= -1e-6);
79		}
80
81	<	/* returns the name of the given tag */
90	<	#ifdef ezxml_name
91	<	#undef ezxml_name
92	<	static char *
93	<	ezxml_name(ezxml_t xml)
94	<	{
95	<	if (xml == NULL)
96	<	return NULL;
97	<	return xml->name;
98	<	}
99	<	#endif
100	<
101	<	/* returns the given tag's character content or empty string if none */
102	<	#ifdef ezxml_txt
103	<	#undef ezxml_txt
104	<	static char *
105	<	ezxml_txt(ezxml_t xml)
106	<	{
107	<	if (xml == NULL)
108	<	return "";
109	<	return xml->txt;
110	<	}
111	<	#endif
112	<
113	<	/* Convert error to standard BSDF code */
81	>	/* convert error to standard BSDF code */
82		static SDError
83		convert_errcode(int ec)
84		{
#	Line 131 | Line 99 | convert_errcode(int ec)
99		return SDEunknown;
100		}
101
102	<	/* Allocate a BSDF matrix of the given size */
102	>	/* allocate a BSDF matrix of the given size */
103		static SDMat *
104		SDnewMatrix(int ni, int no)
105		{
#	Line 155 | Line 123 | SDnewMatrix(int ni, int no)
123		}
124
125		/* Free a BSDF matrix */
126	<	#define SDfreeMatrix            free
126	>	void
127	>	SDfreeMatrix(void *ptr)
128	>	{
129	>	SDMat   *mp = (SDMat *)ptr;
130
131	<	/* get vector for this angle basis index */
132	<	static int
133	<	ab_getvec(FVECT v, int ndx, double randX, void *p)
131	>	if (mp->chroma != NULL) free(mp->chroma);
132	>	free(ptr);
133	>	}
134	>
135	>	/* compute square of real value */
136	>	static double sq(double x) { return x*x; }
137	>
138	>	/* Get vector for this angle basis index (front exiting) */
139	>	int
140	>	fo_getvec(FVECT v, double ndxr, void *p)
141		{
142	<	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
142	>	ANGLE_BASIS     *ab = (ANGLE_BASIS *)p;
143	>	int             ndx = (int)ndxr;
144	>	double          randX = ndxr - ndx;
145		double          rx[2];
146		int             li;
147	<	double          pol, azi, d;
147	>	double          azi, d;
148
149	<	if ((ndx < 0) | (ndx >= ab->nangles))
149	>	if ((ndxr < 0) | (ndx >= ab->nangles))
150		return RC_FAIL;
151		for (li = 0; ndx >= ab->lat[li].nphis; li++)
152		ndx -= ab->lat[li].nphis;
153		SDmultiSamp(rx, 2, randX);
154	<	pol = M_PI/180.*( (1.-rx[0])*ab->lat[li].tmin +
155	<	rx[0]*ab->lat[li+1].tmin );
154	>	d = (1. - rx[0])*sq(cos(M_PI/180.*ab->lat[li].tmin)) +
155	>	rx[0]*sq(cos(M_PI/180.*ab->lat[li+1].tmin));
156	>	v[2] = d = sqrt(d);     /* cos(pol) */
157		azi = 2.*M_PI*(ndx + rx[1] - .5)/ab->lat[li].nphis;
177	–	v[2] = d = cos(pol);
158		d = sqrt(1. - d*d);     /* sin(pol) */
159		v[0] = cos(azi)*d;
160		v[1] = sin(azi)*d;
161		return RC_GOOD;
162		}
163
164	<	/* get index corresponding to the given vector */
165	<	static int
166	<	ab_getndx(const FVECT v, void *p)
164	>	/* Get index corresponding to the given vector (front exiting) */
165	>	int
166	>	fo_getndx(const FVECT v, void *p)
167		{
168	<	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
168	>	ANGLE_BASIS     *ab = (ANGLE_BASIS *)p;
169		int     li, ndx;
170	<	double  pol, azi, d;
170	>	double  pol, azi;
171
172		if (v == NULL)
173		return -1;
174	<	if ((v[2] < .0) | (v[2] > 1.0))
174	>	if ((v[2] < 0) | (v[2] > 1.))
175		return -1;
176	<	pol = 180.0/M_PI*acos(v[2]);
176	>	pol = 180.0/M_PI*Acos(v[2]);
177		azi = 180.0/M_PI*atan2(v[1], v[0]);
178		if (azi < 0.0) azi += 360.0;
179		for (li = 1; ab->lat[li].tmin <= pol; li++)
#	Line 207 | Line 187 | ab_getndx(const FVECT v, void *p)
187		return ndx;
188		}
189
190	<	/* compute square of real value */
191	<	static double sq(double x) { return x*x; }
192	<
213	<	/* get projected solid angle for this angle basis index */
214	<	static double
215	<	ab_getohm(int ndx, void *p)
190	>	/* Get projected solid angle for this angle basis index (universal) */
191	>	double
192	>	io_getohm(int ndx, void *p)
193		{
194	+	static void     *last_p = NULL;
195		static int      last_li = -1;
196		static double   last_ohm;
197		ANGLE_BASIS     *ab = (ANGLE_BASIS *)p;
#	Line 224 | Line 202 | ab_getohm(int ndx, void *p)
202		return -1.;
203		for (li = 0; ndx >= ab->lat[li].nphis; li++)
204		ndx -= ab->lat[li].nphis;
205	<	if (li == last_li)                      /* cached latitude? */
205	>	if ((p == last_p) & (li == last_li))            /* cached latitude? */
206		return last_ohm;
207	+	last_p = p;
208		last_li = li;
230	–	theta1 = M_PI/180. * ab->lat[li+1].tmin;
231	–	if (ab->lat[li].nphis == 1)             /* special case */
232	–	return last_ohm = M_PI*(1. - sq(cos(theta1)));
209		theta = M_PI/180. * ab->lat[li].tmin;
210	+	theta1 = M_PI/180. * ab->lat[li+1].tmin;
211		return last_ohm = M_PI*(sq(cos(theta)) - sq(cos(theta1))) /
212		(double)ab->lat[li].nphis;
213		}
214
215	<	/* get reverse vector for this angle basis index */
216	<	static int
217	<	ab_getvecR(FVECT v, int ndx, double randX, void *p)
215	>	/* Get vector for this angle basis index (back incident) */
216	>	int
217	>	bi_getvec(FVECT v, double ndxr, void *p)
218		{
219	<	int     na = (*(ANGLE_BASIS *)p).nangles;
243	<
244	<	if (!ab_getvec(v, ndx, randX, p))
219	>	if (!fo_getvec(v, ndxr, p))
220		return RC_FAIL;
221
222		v[0] = -v[0];
#	Line 251 | Line 226 | ab_getvecR(FVECT v, int ndx, double randX, void *p)
226		return RC_GOOD;
227		}
228
229	<	/* get index corresponding to the reverse vector */
230	<	static int
231	<	ab_getndxR(const FVECT v, void *p)
229	>	/* Get index corresponding to the vector (back incident) */
230	>	int
231	>	bi_getndx(const FVECT v, void *p)
232		{
233		FVECT  v2;
234
#	Line 261 | Line 236 | ab_getndxR(const FVECT v, void *p)
236		v2[1] = -v[1];
237		v2[2] = -v[2];
238
239	<	return ab_getndx(v2, p);
239	>	return fo_getndx(v2, p);
240		}
241
242	+	/* Get vector for this angle basis index (back exiting) */
243	+	int
244	+	bo_getvec(FVECT v, double ndxr, void *p)
245	+	{
246	+	if (!fo_getvec(v, ndxr, p))
247	+	return RC_FAIL;
248	+
249	+	v[2] = -v[2];
250	+
251	+	return RC_GOOD;
252	+	}
253	+
254	+	/* Get index corresponding to the vector (back exiting) */
255	+	int
256	+	bo_getndx(const FVECT v, void *p)
257	+	{
258	+	FVECT  v2;
259	+
260	+	v2[0] = v[0];
261	+	v2[1] = v[1];
262	+	v2[2] = -v[2];
263	+
264	+	return fo_getndx(v2, p);
265	+	}
266	+
267	+	/* Get vector for this angle basis index (front incident) */
268	+	int
269	+	fi_getvec(FVECT v, double ndxr, void *p)
270	+	{
271	+	if (!fo_getvec(v, ndxr, p))
272	+	return RC_FAIL;
273	+
274	+	v[0] = -v[0];
275	+	v[1] = -v[1];
276	+
277	+	return RC_GOOD;
278	+	}
279	+
280	+	/* Get index corresponding to the vector (front incident) */
281	+	int
282	+	fi_getndx(const FVECT v, void *p)
283	+	{
284	+	FVECT  v2;
285	+
286	+	v2[0] = -v[0];
287	+	v2[1] = -v[1];
288	+	v2[2] = v[2];
289	+
290	+	return fo_getndx(v2, p);
291	+	}
292	+
293	+	/* Get color or grayscale value for BSDF for the given direction pair */
294	+	int
295	+	mBSDF_color(float coef[], const SDMat *dp, int i, int o)
296	+	{
297	+	C_COLOR cxy;
298	+
299	+	coef[0] = mBSDF_value(dp, i, o);
300	+	if (dp->chroma == NULL)
301	+	return 1;       /* grayscale */
302	+
303	+	c_decodeChroma(&cxy, mBSDF_chroma(dp,i,o));
304	+	c_toSharpRGB(&cxy, coef[0], coef);
305	+	coef[0] *= mtx_RGB_coef[0];
306	+	coef[1] *= mtx_RGB_coef[1];
307	+	coef[2] *= mtx_RGB_coef[2];
308	+	return 3;               /* RGB color */
309	+	}
310	+
311		/* load custom BSDF angle basis */
312		static int
313		load_angle_basis(ezxml_t wab)
#	Line 295 | Line 339 | load_angle_basis(ezxml_t wab)
339		ezxml_child(ezxml_child(wbb,
340		"ThetaBounds"), "UpperTheta")));
341		if (!i)
342	<	abase_list[nabases].lat[i].tmin =
299	<	-abase_list[nabases].lat[i+1].tmin;
342	>	abase_list[nabases].lat[0].tmin = 0;
343		else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
344		"ThetaBounds"), "LowerTheta"))),
345		abase_list[nabases].lat[i].tmin)) {
346		sprintf(SDerrorDetail, "Theta values disagree in '%s'",
347	<	abname);
347	>	abname);
348		return RC_DATERR;
349		}
350		abase_list[nabases].nangles +=
#	Line 311 | Line 354 | load_angle_basis(ezxml_t wab)
354		(abase_list[nabases].lat[i].nphis == 1 &&
355		abase_list[nabases].lat[i].tmin > FTINY)) {
356		sprintf(SDerrorDetail, "Illegal phi count in '%s'",
357	<	abname);
357	>	abname);
358		return RC_DATERR;
359		}
360		}
#	Line 357 | Line 400 | get_extrema(SDSpectralDF *df)
400
401		/* load BSDF distribution for this wavelength */
402		static int
403	<	load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
403	>	load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int rowinc)
404		{
405		SDSpectralDF    *df;
406		SDMat           *dp;
#	Line 366 | Line 409 | load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
409		int             i;
410		/* allocate BSDF component */
411		sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
412	<	if (!strcasecmp(sdata, "Transmission Back") || (sd->tf == NULL &&
413	<	!strcasecmp(sdata, "Transmission Front"))) {
414	<	if (sd->tf != NULL)
415	<	SDfreeSpectralDF(sd->tf);
416	<	if ((sd->tf = SDnewSpectralDF(1)) == NULL)
412	>	if (!sdata)
413	>	return RC_FAIL;
414	>	/*
415	>	* Remember that front and back are reversed from WINDOW 6 orientations
416	>	*/
417	>	if (!strcasecmp(sdata, "Transmission Front")) {
418	>	if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(3)) == NULL)
419		return RC_MEMERR;
420	+	df = sd->tb;
421	+	} else if (!strcasecmp(sdata, "Transmission Back")) {
422	+	if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(3)) == NULL)
423	+	return RC_MEMERR;
424		df = sd->tf;
425		} else if (!strcasecmp(sdata, "Reflection Front")) {
426	<	if (sd->rb != NULL)     /* note back-front reversal */
378	<	SDfreeSpectralDF(sd->rb);
379	<	if ((sd->rb = SDnewSpectralDF(1)) == NULL)
426	>	if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(3)) == NULL)
427		return RC_MEMERR;
428		df = sd->rb;
429		} else if (!strcasecmp(sdata, "Reflection Back")) {
430	<	if (sd->rf != NULL)     /* note front-back reversal */
384	<	SDfreeSpectralDF(sd->rf);
385	<	if ((sd->rf = SDnewSpectralDF(1)) == NULL)
430	>	if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(3)) == NULL)
431		return RC_MEMERR;
432		df = sd->rf;
433		} else
434		return RC_FAIL;
435	<	/* XXX should also check "ScatteringDataType" for consistency? */
435	>	/* free previous matrix if any */
436	>	if (df->comp[ct].dist != NULL) {
437	>	SDfreeMatrix(df->comp[ct].dist);
438	>	df->comp[ct].dist = NULL;
439	>	}
440		/* get angle bases */
441		sdata = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
442		if (!sdata || !*sdata) {
#	Line 422 | Line 471 | load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
471		dp->ib_priv = &abase_list[inbi];
472		dp->ob_priv = &abase_list[outbi];
473		if (df == sd->tf) {
474	<	dp->ib_vec = &ab_getvecR;
475	<	dp->ib_ndx = &ab_getndxR;
476	<	dp->ob_vec = &ab_getvec;
477	<	dp->ob_ndx = &ab_getndx;
474	>	dp->ib_vec = &fi_getvec;
475	>	dp->ib_ndx = &fi_getndx;
476	>	dp->ob_vec = &bo_getvec;
477	>	dp->ob_ndx = &bo_getndx;
478	>	} else if (df == sd->tb) {
479	>	dp->ib_vec = &bi_getvec;
480	>	dp->ib_ndx = &bi_getndx;
481	>	dp->ob_vec = &fo_getvec;
482	>	dp->ob_ndx = &fo_getndx;
483		} else if (df == sd->rf) {
484	<	dp->ib_vec = &ab_getvec;
485	<	dp->ib_ndx = &ab_getndx;
486	<	dp->ob_vec = &ab_getvec;
487	<	dp->ob_ndx = &ab_getndx;
484	>	dp->ib_vec = &fi_getvec;
485	>	dp->ib_ndx = &fi_getndx;
486	>	dp->ob_vec = &fo_getvec;
487	>	dp->ob_ndx = &fo_getndx;
488		} else /* df == sd->rb */ {
489	<	dp->ib_vec = &ab_getvecR;
490	<	dp->ib_ndx = &ab_getndxR;
491	<	dp->ob_vec = &ab_getvecR;
492	<	dp->ob_ndx = &ab_getndxR;
489	>	dp->ib_vec = &bi_getvec;
490	>	dp->ib_ndx = &bi_getndx;
491	>	dp->ob_vec = &bo_getvec;
492	>	dp->ob_ndx = &bo_getndx;
493		}
494	<	dp->ib_ohm = &ab_getohm;
495	<	dp->ob_ohm = &ab_getohm;
496	<	df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
497	<	df->comp[0].dist = dp;
444	<	df->comp[0].func = &SDhandleMtx;
494	>	dp->ib_ohm = &io_getohm;
495	>	dp->ob_ohm = &io_getohm;
496	>	df->comp[ct].dist = dp;
497	>	df->comp[ct].func = &SDhandleMtx;
498		/* read BSDF data */
499	<	sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
499	>	sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
500		if (!sdata || !*sdata) {
501		sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
502		sd->name);
503		return RC_FORMERR;
504		}
505		for (i = 0; i < dp->ninc*dp->nout; i++) {
506	<	char  *sdnext = fskip(sdata);
506	>	char    *sdnext = fskip(sdata);
507	>	double  val;
508		if (sdnext == NULL) {
509		sprintf(SDerrorDetail,
510		"Bad/missing BSDF ScatteringData in '%s'",
511		sd->name);
512		return RC_FORMERR;
513		}
514	<	while (*sdnext && isspace(*sdnext))
514	>	while (isspace(*sdnext))
515		sdnext++;
516		if (*sdnext == ',') sdnext++;
517	+	if ((val = atof(sdata)) < 0)
518	+	val = 0;        /* don't allow negative values */
519		if (rowinc) {
520		int     r = i/dp->nout;
521	<	int     c = i - c*dp->nout;
522	<	mBSDF_value(dp,r,c) = atof(sdata);
521	>	int     c = i - r*dp->nout;
522	>	mBSDF_value(dp,r,c) = val;
523		} else
524	<	dp->bsdf[i] = atof(sdata);
524	>	dp->bsdf[i] = val;
525		sdata = sdnext;
526		}
527	<	return get_extrema(df);
527	>	return (ct == mtx_Y) ? get_extrema(df) : RC_GOOD;
528		}
529
530	<	/* Subtract minimum (diffuse) scattering amount from BSDF */
530	>	/* copy our RGB (x,y) primary chromaticities */
531	>	static void
532	>	copy_RGB_prims(C_COLOR cspec[])
533	>	{
534	>	if (mtx_RGB_coef[1] < .001) {   /* need to initialize */
535	>	int     i = 3;
536	>	while (i--) {
537	>	float   rgb[3];
538	>	rgb[0] = rgb[1] = rgb[2] = .0f;
539	>	rgb[i] = 1.f;
540	>	mtx_RGB_coef[i] = c_fromSharpRGB(rgb, &mtx_RGB_prim[i]);
541	>	}
542	>	}
543	>	memcpy(cspec, mtx_RGB_prim, sizeof(mtx_RGB_prim));
544	>	}
545	>
546	>	/* encode chromaticity if XYZ -- reduce to one channel in any case */
547	>	static SDSpectralDF *
548	>	encode_chroma(SDSpectralDF *df)
549	>	{
550	>	SDMat   *mpx, *mpy, *mpz;
551	>	int     n;
552	>
553	>	if (df == NULL || df->ncomp != 3)
554	>	return df;
555	>
556	>	mpy = (SDMat *)df->comp[mtx_Y].dist;
557	>	if (mpy == NULL) {
558	>	free(df);
559	>	return NULL;
560	>	}
561	>	mpx = (SDMat *)df->comp[mtx_X].dist;
562	>	mpz = (SDMat *)df->comp[mtx_Z].dist;
563	>	if (mpx == NULL || (mpx->ninc != mpy->ninc) | (mpx->nout != mpy->nout))
564	>	goto done;
565	>	if (mpz == NULL || (mpz->ninc != mpy->ninc) | (mpz->nout != mpy->nout))
566	>	goto done;
567	>	mpy->chroma = (C_CHROMA *)malloc(sizeof(C_CHROMA)*mpy->ninc*mpy->nout);
568	>	if (mpy->chroma == NULL)
569	>	goto done;              /* XXX punt */
570	>	/* encode chroma values */
571	>	for (n = mpy->ninc*mpy->nout; n--; ) {
572	>	const double    sum = mpx->bsdf[n] + mpy->bsdf[n] + mpz->bsdf[n];
573	>	C_COLOR         cxy;
574	>	if (sum > .0)
575	>	c_cset(&cxy, mpx->bsdf[n]/sum, mpy->bsdf[n]/sum);
576	>	else
577	>	c_cset(&cxy, 1./3., 1./3.);
578	>	mpy->chroma[n] = c_encodeChroma(&cxy);
579	>	}
580	>	done:                                   /* free X & Z channels */
581	>	if (mpx != NULL) SDfreeMatrix(mpx);
582	>	if (mpz != NULL) SDfreeMatrix(mpz);
583	>	if (mpy->chroma == NULL)        /* grayscale after all? */
584	>	df->comp[0].cspec[0] = c_dfcolor;
585	>	else                            /* else copy RGB primaries */
586	>	copy_RGB_prims(df->comp[0].cspec);
587	>	df->ncomp = 1;                  /* return resized struct */
588	>	return (SDSpectralDF *)realloc(df, sizeof(SDSpectralDF));
589	>	}
590	>
591	>	/* subtract minimum (diffuse) scattering amount from BSDF */
592		static double
593	<	subtract_min(SDMat *sm)
593	>	subtract_min(C_COLOR *cs, SDMat *sm)
594		{
595	<	float   minv = sm->bsdf[0];
596	<	int     n = sm->ninc*sm->nout;
597	<	int     i;
595	>	const int       ncomp = 1 + 2*(sm->chroma != NULL);
596	>	float           min_coef[3], ymin, coef[3];
597	>	int             i, o, c;
598
599	<	for (i = n; --i; )
600	<	if (sm->bsdf[i] < minv)
601	<	minv = sm->bsdf[i];
602	<	for (i = n; i--; )
603	<	sm->bsdf[i] -= minv;
599	>	min_coef[0] = min_coef[1] = min_coef[2] = FHUGE;
600	>	for (i = 0; i < sm->ninc; i++)
601	>	for (o = 0; o < sm->nout; o++) {
602	>	c = mBSDF_color(coef, sm, i, o);
603	>	while (c--)
604	>	if (coef[c] < min_coef[c])
605	>	min_coef[c] = coef[c];
606	>	}
607	>	ymin = 0;
608	>	for (c = ncomp; c--; )
609	>	ymin += min_coef[c];
610	>	if (ymin <= .01/M_PI)           /* not worth bothering about? */
611	>	return .0;
612	>	if (ncomp == 1) {               /* subtract grayscale minimum */
613	>	for (i = sm->ninc*sm->nout; i--; )
614	>	sm->bsdf[i] -= ymin;
615	>	*cs = c_dfcolor;
616	>	return M_PI*ymin;
617	>	}
618	>	/* else subtract colored minimum */
619	>	for (i = 0; i < sm->ninc; i++)
620	>	for (o = 0; o < sm->nout; o++) {
621	>	C_COLOR cxy;
622	>	c = mBSDF_color(coef, sm, i, o);
623	>	while (c--)
624	>	coef[c] = (coef[c] - min_coef[c]) /
625	>	mtx_RGB_coef[c];
626	>	if (c_fromSharpRGB(coef, &cxy) > 1e-5)
627	>	mBSDF_chroma(sm,i,o) = c_encodeChroma(&cxy);
628	>	mBSDF_value(sm,i,o) -= ymin;
629	>	}
630	>	/* return colored minimum */
631	>	for (i = 3; i--; )
632	>	coef[i] = min_coef[i]/mtx_RGB_coef[i];
633	>	c_fromSharpRGB(coef, cs);
634
635	<	return minv*M_PI;               /* be sure to include multiplier */
635	>	return M_PI*ymin;
636		}
637
638	<	/* Extract and separate diffuse portion of BSDF */
639	<	static void
638	>	/* Extract and separate diffuse portion of BSDF & convert color */
639	>	static SDSpectralDF *
640		extract_diffuse(SDValue *dv, SDSpectralDF *df)
641		{
495	–	int     n;
642
643	+	df = encode_chroma(df);         /* reduce XYZ to Y + chroma */
644		if (df == NULL || df->ncomp <= 0) {
645		dv->spec = c_dfcolor;
646		dv->cieY = .0;
647	<	return;
647	>	return df;
648		}
649	<	dv->spec = df->comp[0].cspec[0];
650	<	dv->cieY = subtract_min((SDMat *)df->comp[0].dist);
651	<	/* in case of multiple components */
505	<	for (n = df->ncomp; --n; ) {
506	<	double  ymin = subtract_min((SDMat *)df->comp[n].dist);
507	<	c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
508	<	dv->cieY += ymin;
509	<	}
510	<	df->maxHemi -= dv->cieY;        /* adjust minimum hemispherical */
649	>	/* subtract minimum value */
650	>	dv->cieY = subtract_min(&dv->spec, (SDMat *)df->comp[0].dist);
651	>	df->maxHemi -= dv->cieY;        /* adjust maximum hemispherical */
652		/* make sure everything is set */
653		c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
654	+	return df;
655		}
656
657		/* Load a BSDF matrix from an open XML file */
658		SDError
659		SDloadMtx(SDData *sd, ezxml_t wtl)
660		{
661	<	ezxml_t                 wld, wdb;
662	<	int                     rowIn;
663	<	struct BSDF_data        *dp;
664	<	char                    *txt;
665	<	int                     rval;
524	<
661	>	ezxml_t         wld, wdb;
662	>	int             rowIn;
663	>	char            *txt;
664	>	int             rval;
665	>	/* basic checks and data ordering */
666		txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
667		"DataDefinition"), "IncidentDataStructure"));
668		if (txt == NULL || !*txt) {
#	Line 540 | Line 681 | SDloadMtx(SDData *sd, ezxml_t wtl)
681		sd->name);
682		return SDEsupport;
683		}
684	<	/* get angle basis */
685	<	rval = load_angle_basis(ezxml_child(ezxml_child(wtl,
686	<	"DataDefinition"), "AngleBasis"));
687	<	if (rval < 0)
688	<	return convert_errcode(rval);
689	<	/* load BSDF components */
684	>	/* get angle bases */
685	>	for (wld = ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis");
686	>	wld != NULL; wld = wld->next) {
687	>	rval = load_angle_basis(wld);
688	>	if (rval < 0)
689	>	return convert_errcode(rval);
690	>	}
691	>	/* load BSDF components */
692		for (wld = ezxml_child(wtl, "WavelengthData");
693		wld != NULL; wld = wld->next) {
694	<	if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
695	<	"Visible"))
696	<	continue;       /* just visible for now */
694	>	const char      *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
695	>	int             ct = -1;
696	>	if (!strcasecmp(cnm, "Visible"))
697	>	ct = mtx_Y;
698	>	else if (!strcasecmp(cnm, "CIE-X"))
699	>	ct = mtx_X;
700	>	else if (!strcasecmp(cnm, "CIE-Z"))
701	>	ct = mtx_Z;
702	>	else
703	>	continue;
704		for (wdb = ezxml_child(wld, "WavelengthDataBlock");
705		wdb != NULL; wdb = wdb->next)
706	<	if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0)
706	>	if ((rval = load_bsdf_data(sd, wdb, ct, rowIn)) < 0)
707		return convert_errcode(rval);
708		}
709	<	/* separate diffuse components */
710	<	extract_diffuse(&sd->rLambFront, sd->rf);
711	<	extract_diffuse(&sd->rLambBack, sd->rb);
712	<	extract_diffuse(&sd->tLamb, sd->tf);
713	<	/* return success */
709	>	/* separate diffuse components */
710	>	sd->rf = extract_diffuse(&sd->rLambFront, sd->rf);
711	>	sd->rb = extract_diffuse(&sd->rLambBack, sd->rb);
712	>	if (sd->tf != NULL)
713	>	sd->tf = extract_diffuse(&sd->tLamb, sd->tf);
714	>	if (sd->tb != NULL)
715	>	sd->tb = extract_diffuse(&sd->tLamb, sd->tb);
716	>	/* return success */
717		return SDEnone;
718		}
719
720		/* Get Matrix BSDF value */
721		static int
722		SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
723	<	const FVECT inVec, const void *dist)
723	>	const FVECT inVec, SDComponent *sdc)
724		{
725	<	const SDMat     *dp = (const SDMat *)dist;
725	>	const SDMat     *dp;
726		int             i_ndx, o_ndx;
727	+	/* check arguments */
728	+	if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
729	+	|| (dp = (SDMat *)sdc->dist) == NULL)
730	+	return 0;
731		/* get angle indices */
732		i_ndx = mBSDF_incndx(dp, inVec);
733		o_ndx = mBSDF_outndx(dp, outVec);
#	Line 581 | Line 738 | SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
738		}
739		if ((i_ndx < 0) | (o_ndx < 0))
740		return 0;               /* nothing from this component */
741	<	coef[0] = mBSDF_value(dp, i_ndx, o_ndx);
742	<	return 1;                       /* XXX monochrome for now */
741	>
742	>	return mBSDF_color(coef, dp, i_ndx, o_ndx);
743		}
744
745	<	/* Query solid angle for vector */
745	>	/* Query solid angle for vector(s) */
746		static SDError
747	<	SDqueryMtxProjSA(double *psa, const FVECT vec, int qflags, const void *dist)
747	>	SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2,
748	>	int qflags, SDComponent *sdc)
749		{
750	<	const SDMat     *dp = (const SDMat *)dist;
750	>	const SDMat     *dp;
751		double          inc_psa, out_psa;
752		/* check arguments */
753	<	if ((psa == NULL) | (vec == NULL) | (dp == NULL))
753	>	if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
754	>	(dp = (SDMat *)sdc->dist) == NULL)
755		return SDEargument;
756	+	if (v2 == NULL)
757	+	v2 = v1;
758		/* get projected solid angles */
759	<	inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec));
760	<	out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec));
759	>	out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1));
760	>	inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2));
761	>	if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) {
762	>	inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2));
763	>	out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1));
764	>	}
765
766		switch (qflags) {               /* record based on flag settings */
602	–	case SDqueryVal:
603	–	psa[0] = .0;
604	–	/* fall through */
767		case SDqueryMax:
768		if (inc_psa > psa[0])
769		psa[0] = inc_psa;
#	Line 609 | Line 771 | SDqueryMtxProjSA(double *psa, const FVECT vec, int qfl
771		psa[0] = out_psa;
772		break;
773		case SDqueryMin+SDqueryMax:
774	<	if (inc_psa > psa[0])
774	>	if (inc_psa > psa[1])
775		psa[1] = inc_psa;
776	<	if (out_psa > psa[0])
776	>	if (out_psa > psa[1])
777		psa[1] = out_psa;
778		/* fall through */
779	+	case SDqueryVal:
780	+	if (qflags == SDqueryVal)
781	+	psa[0] = M_PI;
782	+	/* fall through */
783		case SDqueryMin:
784	<	if ((inc_psa > .0) & (inc_psa < psa[0]))
784	>	if ((inc_psa > 0) & (inc_psa < psa[0]))
785		psa[0] = inc_psa;
786	<	if ((out_psa > .0) & (out_psa < psa[0]))
786	>	if ((out_psa > 0) & (out_psa < psa[0]))
787		psa[0] = out_psa;
788		break;
789		}
790		/* make sure it's legal */
791	<	return (psa[0] <= .0) ? SDEinternal : SDEnone;
791	>	return (psa[0] <= 0) ? SDEinternal : SDEnone;
792		}
793
794		/* Compute new cumulative distribution from BSDF */
#	Line 660 | Line 826 | make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp
826		static const SDCDst *
827		SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
828		{
829	<	SDMat           *dp = (SDMat *)sdc->dist;
829	>	SDMat           *dp;
830		int             reverse;
831		SDMatCDst       myCD;
832		SDMatCDst       *cd, *cdlast;
833		/* check arguments */
834	<	if ((inVec == NULL) | (dp == NULL))
834	>	if ((inVec == NULL) | (sdc == NULL) ||
835	>	(dp = (SDMat *)sdc->dist) == NULL)
836		return NULL;
837		memset(&myCD, 0, sizeof(myCD));
838		myCD.indx = mBSDF_incndx(dp, inVec);
#	Line 684 | Line 851 | SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
851		reverse = 1;
852		}
853		cdlast = NULL;                  /* check for it in cache list */
854	<	for (cd = (SDMatCDst *)sdc->cdList;
855	<	cd != NULL; cd = (SDMatCDst *)cd->next) {
854	>	/* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
855	>	for (cd = (SDMatCDst *)sdc->cdList; cd != NULL;
856	>	cdlast = cd, cd = cd->next)
857		if (cd->indx == myCD.indx && (cd->calen == myCD.calen) &
858		(cd->ob_priv == myCD.ob_priv) &
859		(cd->ob_vec == myCD.ob_vec))
860		break;
693	–	cdlast = cd;
694	–	}
861		if (cd == NULL) {               /* need to allocate new entry */
862		cd = (SDMatCDst *)malloc(sizeof(SDMatCDst) +
863	<	myCD.calen*sizeof(myCD.carr[0]));
863	>	sizeof(myCD.carr[0])*myCD.calen);
864		if (cd == NULL)
865		return NULL;
866		*cd = myCD;             /* compute cumulative distribution */
#	Line 706 | Line 872 | SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
872		}
873		if (cdlast != NULL) {           /* move entry to head of cache list */
874		cdlast->next = cd->next;
875	<	cd->next = sdc->cdList;
875	>	cd->next = (SDMatCDst *)sdc->cdList;
876		sdc->cdList = (SDCDst *)cd;
877		}
878	+	/* END MUTEX LOCK */
879		return (SDCDst *)cd;            /* ready to go */
880		}
881
882		/* Sample cumulative distribution */
883		static SDError
884	<	SDsampMtxCDist(FVECT outVec, double randX, const SDCDst *cdp)
884	>	SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp)
885		{
886		const unsigned  maxval = ~0;
887		const SDMatCDst *mcd = (const SDMatCDst *)cdp;
888		const unsigned  target = randX*maxval;
889		int             i, iupper, ilower;
890		/* check arguments */
891	<	if ((outVec == NULL) | (mcd == NULL))
891	>	if ((ioVec == NULL) | (mcd == NULL))
892		return SDEargument;
893		/* binary search to find index */
894		ilower = 0; iupper = mcd->calen;
895		while ((i = (iupper + ilower) >> 1) != ilower)
896	<	if ((long)target >= (long)mcd->carr[i])
896	>	if (target >= mcd->carr[i])
897		ilower = i;
898		else
899		iupper = i;
#	Line 734 | Line 901 | SDsampMtxCDist(FVECT outVec, double randX, const SDCDs
901		randX = (randX*maxval - mcd->carr[ilower]) /
902		(double)(mcd->carr[iupper] - mcd->carr[ilower]);
903		/* convert index to vector */
904	<	if ((*mcd->ob_vec)(outVec, i, randX, mcd->ob_priv))
904	>	if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv))
905		return SDEnone;
906	<	strcpy(SDerrorDetail, "BSDF sampling fault");
906	>	strcpy(SDerrorDetail, "Matrix BSDF sampling fault");
907		return SDEinternal;
908		}
909
910		/* Fixed resolution BSDF methods */
911	<	SDFunc                  SDhandleMtx = {
911	>	const SDFunc            SDhandleMtx = {
912		&SDgetMtxBSDF,
913		&SDqueryMtxProjSA,
914		&SDgetMtxCDist,

Diff Legend

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