ViewVC Help

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

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

Comparing ray/src/common/bsdf_t.c (file contents):
Revision 3.30 by greg, Thu Jul 4 15:14:45 2013 UTC vs.
Revision 3.40 by greg, Mon Apr 6 18:28:35 2015 UTC

#	Line 21 | Line 21 | static const char RCSid[] = "$Id$";
21		#include "hilbert.h"
22
23		/* Callback function type for SDtraverseTre() */
24	<	typedef int     SDtreCallback(float val, const double *cmin,
25	<	double csiz, void *cptr);
26	<
24	>	typedef int     SDtreCallback(float val, const double *cmin, double csiz,
25	>	void *cptr);
26		/* reference width maximum (1.0) */
27		static const unsigned   iwbits = sizeof(unsigned)*4;
28		static const unsigned   iwmax = 1<<(sizeof(unsigned)*4);
#	Line 33 | Line 32 | static const unsigned  cumlmax = ~0;
32		static const FVECT      zvec = {.0, .0, 1.};
33		/* quantization value */
34		static double           quantum = 1./256.;
35	+	/* our RGB primaries */
36	+	static C_COLOR          tt_RGB_prim[3];
37	+	static float            tt_RGB_coef[3];
38
39	+	static const double     czero[SD_MAXDIM];
40	+
41	+	enum {tt_Y, tt_u, tt_v};                /* tree components (tt_Y==0) */
42	+
43		/* Struct used for our distribution-building callback */
44		typedef struct {
45		short           nic;            /* number of input coordinates */
#	Line 108 | Line 114 | SDFreeBTre(void *p)
114
115		if (sdt == NULL)
116		return;
117	<	SDfreeTre(sdt->st);
117	>	SDfreeTre(sdt->stc[tt_Y]);
118	>	SDfreeTre(sdt->stc[tt_u]);
119	>	SDfreeTre(sdt->stc[tt_v]);
120		free(sdt);
121		}
122
#	Line 171 | Line 179 | SDsimplifyTre(SDNode *st)
179		return st;
180		}
181
182	+	/* Assign the given voxel in tree (produces no grid nodes) */
183	+	static SDNode *
184	+	SDsetVoxel(SDNode *sroot, int nd, const double *tmin, const double tsiz, float val)
185	+	{
186	+	double  ctrk[SD_MAXDIM];
187	+	double  csiz = 1.;
188	+	SDNode  *st;
189	+	int     i, n;
190	+	/* check arguments */
191	+	for (i = nd; i-- > 0; )
192	+	if ((tmin[i] < .0) | (tmin[i] >= 1.-FTINY))
193	+	break;
194	+	if ((i >= 0) | (nd <= 0) | (tsiz <= FTINY) | (tsiz > 1.+FTINY) |
195	+	(sroot != NULL && sroot->ndim != nd)) {
196	+	SDfreeTre(sroot);
197	+	return NULL;
198	+	}
199	+	if (tsiz >= 1.-FTINY) {         /* special case when tree is a leaf */
200	+	SDfreeTre(sroot);
201	+	if ((sroot = SDnewNode(nd, 0)) != NULL)
202	+	sroot->u.v[0] = val;
203	+	return sroot;
204	+	}
205	+	/* make sure we have branching root */
206	+	if (sroot != NULL && sroot->log2GR >= 0) {
207	+	SDfreeTre(sroot); sroot = NULL;
208	+	}
209	+	if (sroot == NULL && (sroot = SDnewNode(nd, -1)) == NULL)
210	+	return NULL;
211	+	st = sroot;                     /* climb/grow tree */
212	+	memset(ctrk, 0, sizeof(ctrk));
213	+	for ( ; ; ) {
214	+	csiz *= .5;             /* find appropriate branch */
215	+	n = 0;
216	+	for (i = nd; i--; )
217	+	if (ctrk[i]+csiz <= tmin[i]+FTINY) {
218	+	ctrk[i] += csiz;
219	+	n |= 1 << i;
220	+	}
221	+	/* reached desired voxel? */
222	+	if (csiz <= tsiz+FTINY) {
223	+	SDfreeTre(st->u.t[n]);
224	+	st = st->u.t[n] = SDnewNode(nd, 0);
225	+	break;
226	+	}
227	+	/* else grow tree as needed */
228	+	if (st->u.t[n] != NULL && st->u.t[n]->log2GR >= 0) {
229	+	SDfreeTre(st->u.t[n]); st->u.t[n] = NULL;
230	+	}
231	+	if (st->u.t[n] == NULL)
232	+	st->u.t[n] = SDnewNode(nd, -1);
233	+	if ((st = st->u.t[n]) == NULL)
234	+	break;
235	+	}
236	+	if (st == NULL) {
237	+	SDfreeTre(sroot);
238	+	return NULL;
239	+	}
240	+	st->u.v[0] = val;               /* assign leaf and return root */
241	+	return sroot;
242	+	}
243	+
244		/* Find smallest leaf in tree */
245		static double
246		SDsmallestLeaf(const SDNode *st)
#	Line 279 | Line 349 | SDdotravTre(const SDNode *st, const double *pos, int c
349		int     rv, rval = 0;
350		double  bmin[SD_MAXDIM];
351		int     i, n;
352	+	/* paranoia */
353	+	if (st == NULL)
354	+	return 0;
355		/* in branches? */
356		if (st->log2GR < 0) {
357		unsigned        skipmask = 0;
358		csiz *= .5;
359		for (i = st->ndim; i--; )
360	<	if (1<<i & cmask)
360	>	if (1<<i & cmask) {
361		if (pos[i] < cmin[i] + csiz)
362		for (n = 1 << st->ndim; n--; ) {
363		if (n & 1<<i)
#	Line 295 | Line 368 | SDdotravTre(const SDNode *st, const double *pos, int c
368		if (!(n & 1<<i))
369		skipmask |= 1<<n;
370		}
371	+	}
372		for (n = 1 << st->ndim; n--; ) {
373		if (1<<n & skipmask)
374		continue;
#	Line 376 | Line 450 | static int
450		SDtraverseTre(const SDNode *st, const double *pos, int cmask,
451		SDtreCallback *cf, void *cptr)
452		{
379	–	static double   czero[SD_MAXDIM];
453		int             i;
454		/* check arguments */
455		if ((st == NULL) | (cf == NULL))
#	Line 401 | Line 474 | SDlookupTre(const SDNode *st, const double *pos, doubl
474		hcube[i] = .0;
475		}
476		/* climb the tree */
477	<	while (st->log2GR < 0) {
477	>	while (st != NULL && st->log2GR < 0) {
478		n = 0;                  /* move to appropriate branch */
479		if (hcube) hcube[st->ndim] *= .5;
480		for (i = st->ndim; i--; ) {
#	Line 414 | Line 487 | SDlookupTre(const SDNode *st, const double *pos, doubl
487		st = st->u.t[n];        /* avoids tail recursion */
488		pos = spos;
489		}
490	+	if (st == NULL)                 /* should never happen? */
491	+	return .0;
492		if (st->log2GR == 0)            /* short cut */
493		return st->u.v[0];
494		n = t = 0;                      /* find grid array index */
#	Line 429 | Line 504 | SDlookupTre(const SDNode *st, const double *pos, doubl
504		return st->u.v[n];              /* no interpolation */
505		}
506
507	+	/* Convert CIE (Y,u',v') color to our RGB */
508	+	static void
509	+	SDyuv2rgb(double yval, double uprime, double vprime, float rgb[3])
510	+	{
511	+	const double    dfact = 1./(6.*uprime - 16.*vprime + 12.);
512	+	C_COLOR         cxy;
513	+
514	+	c_cset(&cxy, 9.*uprime*dfact, 4.*vprime*dfact);
515	+	c_toSharpRGB(&cxy, yval, rgb);
516	+	}
517	+
518		/* Query BSDF value and sample hypercube for the given vectors */
519	<	static float
520	<	SDqueryTre(const SDTre *sdt, const FVECT outVec, const FVECT inVec, double *hc)
519	>	static int
520	>	SDqueryTre(const SDTre *sdt, float *coef,
521	>	const FVECT outVec, const FVECT inVec, double *hc)
522		{
523		const RREAL     *vtmp;
524	+	float           yval;
525		FVECT           rOutVec;
526		double          gridPos[4];
527
528	+	if (sdt->stc[tt_Y] == NULL)     /* paranoia, I hope */
529	+	return 0;
530	+
531		switch (sdt->sidef) {           /* whose side are you on? */
532		case SD_FREFL:
533		if ((outVec[2] < 0) | (inVec[2] < 0))
534	<	return -1.;
534	>	return 0;
535		break;
536		case SD_BREFL:
537		if ((outVec[2] > 0) | (inVec[2] > 0))
538	<	return -1.;
538	>	return 0;
539		break;
540		case SD_FXMIT:
541		if (outVec[2] > 0) {
542		if (inVec[2] > 0)
543	<	return -1.;
543	>	return 0;
544		vtmp = outVec; outVec = inVec; inVec = vtmp;
545		} else if (inVec[2] < 0)
546	<	return -1.;
546	>	return 0;
547		break;
548		case SD_BXMIT:
549		if (inVec[2] > 0) {
550		if (outVec[2] > 0)
551	<	return -1.;
551	>	return 0;
552		vtmp = outVec; outVec = inVec; inVec = vtmp;
553		} else if (outVec[2] < 0)
554	<	return -1.;
554	>	return 0;
555		break;
556		default:
557	<	return -1.;
557	>	return 0;
558		}
559		/* convert vector coordinates */
560	<	if (sdt->st->ndim == 3) {
560	>	if (sdt->stc[tt_Y]->ndim == 3) {
561		spinvector(rOutVec, outVec, zvec, -atan2(-inVec[1],-inVec[0]));
562	<	gridPos[0] = .5 - .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
562	>	gridPos[0] = (.5-FTINY) -
563	>	.5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
564		SDdisk2square(gridPos+1, rOutVec[0], rOutVec[1]);
565	<	} else if (sdt->st->ndim == 4) {
565	>	} else if (sdt->stc[tt_Y]->ndim == 4) {
566		SDdisk2square(gridPos, -inVec[0], -inVec[1]);
567		SDdisk2square(gridPos+2, outVec[0], outVec[1]);
568		} else
569	<	return -1.;             /* should be internal error */
570	<
571	<	return SDlookupTre(sdt->st, gridPos, hc);
569	>	return 0;               /* should be internal error */
570	>	/* get BSDF value */
571	>	yval = SDlookupTre(sdt->stc[tt_Y], gridPos, hc);
572	>	if (sdt->stc[tt_u] == NULL || sdt->stc[tt_v] == NULL) {
573	>	if (coef != NULL) *coef = yval;
574	>	return 1;               /* no color */
575	>	}
576	>	if (coef == NULL)               /* just getting hypercube? */
577	>	return 1;
578	>	/* else decode color */
579	>	SDyuv2rgb(yval, SDlookupTre(sdt->stc[tt_u], gridPos, NULL),
580	>	SDlookupTre(sdt->stc[tt_v], gridPos, NULL), coef);
581	>	coef[0] *= tt_RGB_coef[0];
582	>	coef[1] *= tt_RGB_coef[1];
583	>	coef[2] *= tt_RGB_coef[2];
584	>	return 3;
585		}
586
587		/* Compute non-diffuse component for variable-resolution BSDF */
#	Line 489 | Line 594 | SDgetTreBSDF(float coef[SDmaxCh], const FVECT outVec,
594		|| sdc->dist == NULL)
595		return 0;
596		/* get nearest BSDF value */
597	<	coef[0] = SDqueryTre((SDTre *)sdc->dist, outVec, inVec, NULL);
493	<	return (coef[0] >= 0);          /* monochromatic for now */
597	>	return SDqueryTre((SDTre *)sdc->dist, coef, outVec, inVec, NULL);
598		}
599
600		/* Callback to build cumulative distribution using SDtraverseTre() */
#	Line 515 | Line 619 | build_scaffold(float val, const double *cmin, double c
619		sp->wmax = wid;
620		if (sp->alen >= sp->nall) {     /* need more space? */
621		struct outdir_s *ndarr;
622	<	sp->nall += 1024;
622	>	sp->nall = (int)(1.5*sp->nall) + 256;
623		ndarr = (struct outdir_s *)realloc(sp->darr,
624		sizeof(struct outdir_s)*sp->nall);
625		if (ndarr == NULL) {
#	Line 566 | Line 670 | make_cdist(const SDTre *sdt, const double *invec, int
670		/* initialize scaffold */
671		myScaffold.wmin = iwmax;
672		myScaffold.wmax = 0;
673	<	myScaffold.nic = sdt->st->ndim - 2;
673	>	myScaffold.nic = sdt->stc[tt_Y]->ndim - 2;
674		myScaffold.rev = rev;
675		myScaffold.alen = 0;
676		myScaffold.nall = 512;
#	Line 580 | Line 684 | make_cdist(const SDTre *sdt, const double *invec, int
684		pos[i+2*rev] = invec[i];
685		cmask <<= 2*rev;
686		/* grow the distribution */
687	<	if (SDtraverseTre(sdt->st, pos, cmask,
688	<	&build_scaffold, &myScaffold) < 0) {
687	>	if (SDtraverseTre(sdt->stc[tt_Y], pos, cmask,
688	>	build_scaffold, &myScaffold) < 0) {
689		free(myScaffold.darr);
690		return NULL;
691		}
#	Line 598 | Line 702 | make_cdist(const SDTre *sdt, const double *invec, int
702		cd->isodist = (myScaffold.nic == 1);
703		/* sort the distribution */
704		qsort(myScaffold.darr, cd->calen = myScaffold.alen,
705	<	sizeof(struct outdir_s), &sscmp);
705	>	sizeof(struct outdir_s), sscmp);
706
707		/* record input range */
708		scale = myScaffold.wmin / (double)iwmax;
#	Line 669 | Line 773 | SDgetTreCDist(const FVECT inVec, SDComponent *sdc)
773		default:
774		return NULL;
775		}
776	<	if (sdt->st->ndim == 3) {       /* isotropic BSDF? */
776	>	if (sdt->stc[tt_Y]->ndim == 3) {        /* isotropic BSDF? */
777		if (mode != sdt->sidef) /* XXX unhandled reciprocity */
778		return &SDemptyCD;
779	<	inCoord[0] = .5 - .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
780	<	} else if (sdt->st->ndim == 4) {
779	>	inCoord[0] = (.5-FTINY) -
780	>	.5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
781	>	} else if (sdt->stc[tt_Y]->ndim == 4) {
782		if (mode != sdt->sidef) /* use reciprocity? */
783		SDdisk2square(inCoord, inVec[0], inVec[1]);
784		else
#	Line 681 | Line 786 | SDgetTreCDist(const FVECT inVec, SDComponent *sdc)
786		} else
787		return NULL;            /* should be internal error */
788		/* quantize to avoid f.p. errors */
789	<	for (i = sdt->st->ndim - 2; i--; )
789	>	for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
790		inCoord[i] = floor(inCoord[i]/quantum)*quantum + .5*quantum;
791		cdlast = NULL;                  /* check for direction in cache list */
792		for (cd = (SDTreCDst *)sdc->cdList; cd != NULL;
793		cdlast = cd, cd = cd->next) {
794		if (cd->sidef != mode)
795		continue;
796	<	for (i = sdt->st->ndim - 2; i--; )
796	>	for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
797		if ((cd->clim[i][0] > inCoord[i]) |
798		(inCoord[i] >= cd->clim[i][1]))
799		break;
#	Line 718 | Line 823 | SDqueryTreProjSA(double *psa, const FVECT v1, const RR
823		/* get projected solid angle(s) */
824		if (v2 != NULL) {
825		const SDTre     *sdt = (SDTre *)sdc->dist;
826	<	double          hcube[SD_MAXDIM];
827	<	if (SDqueryTre(sdt, v1, v2, hcube) < 0) {
826	>	double          hcube[SD_MAXDIM+1];
827	>	if (!SDqueryTre(sdt, NULL, v1, v2, hcube)) {
828		strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA");
829		return SDEinternal;
830		}
831	<	myPSA[0] = hcube[sdt->st->ndim];
831	>	myPSA[0] = hcube[sdt->stc[tt_Y]->ndim];
832		myPSA[1] = myPSA[0] *= myPSA[0] * M_PI;
833		} else {
834		const SDTreCDst *cd = (const SDTreCDst *)SDgetTreCDist(v1, sdc);
#	Line 797 | Line 902 | SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
902		SDsquare2disk(gpos, gpos[0], gpos[1]);
903		/* compute Z-coordinate */
904		gpos[2] = 1. - gpos[0]*gpos[0] - gpos[1]*gpos[1];
905	<	if (gpos[2] > 0)                /* paranoia, I hope */
801	<	gpos[2] = sqrt(gpos[2]);
905	>	gpos[2] = sqrt(gpos[2]*(gpos[2]>0));
906		/* emit from back? */
907		if ((cd->sidef == SD_BREFL) | (cd->sidef == SD_FXMIT))
908		gpos[2] = -gpos[2];
909	<	if (cd->isodist) {              /* rotate isotropic result */
909	>	if (cd->isodist) {              /* rotate isotropic sample */
910		rotangle = atan2(-ioVec[1],-ioVec[0]);
911	<	VCOPY(ioVec, gpos);
808	<	spinvector(ioVec, ioVec, zvec, rotangle);
911	>	spinvector(ioVec, gpos, zvec, rotangle);
912		} else
913		VCOPY(ioVec, gpos);
914		return SDEnone;
#	Line 847 | Line 950 | load_values(char **spp, float *va, int n)
950		char    *svnext;
951
952		while (n-- > 0 && (svnext = fskip(*spp)) != NULL) {
953	<	*v++ = atof(*spp);
953	>	if ((*v++ = atof(*spp)) < 0)
954	>	v[-1] = 0;
955		*spp = svnext;
956		eat_token(spp, ',');
957		}
#	Line 906 | Line 1010 | load_tree_data(char **spp, int nd)
1010		static SDError
1011		get_extrema(SDSpectralDF *df)
1012		{
1013	<	SDNode  *st = (*(SDTre *)df->comp[0].dist).st;
1013	>	SDNode  *st = (*(SDTre *)df->comp[0].dist).stc[tt_Y];
1014		double  stepWidth, dhemi, bmin[4], bmax[4];
1015
1016		stepWidth = SDsmallestLeaf(st);
#	Line 946 | Line 1050 | get_extrema(SDSpectralDF *df)
1050
1051		/* Load BSDF distribution for this wavelength */
1052		static SDError
1053	<	load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
1053	>	load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int ndim)
1054		{
1055		SDSpectralDF    *df;
1056		SDTre           *sdt;
#	Line 959 | Line 1063 | load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
1063		* Remember that front and back are reversed from WINDOW 6 orientations
1064		*/
1065		if (!strcasecmp(sdata, "Transmission Front")) {
1066	<	if (sd->tb != NULL)
963	<	SDfreeSpectralDF(sd->tb);
964	<	if ((sd->tb = SDnewSpectralDF(1)) == NULL)
1066	>	if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(1)) == NULL)
1067		return SDEmemory;
1068		df = sd->tb;
1069		} else if (!strcasecmp(sdata, "Transmission Back")) {
1070	<	if (sd->tf != NULL)
969	<	SDfreeSpectralDF(sd->tf);
970	<	if ((sd->tf = SDnewSpectralDF(1)) == NULL)
1070	>	if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(1)) == NULL)
1071		return SDEmemory;
1072		df = sd->tf;
1073		} else if (!strcasecmp(sdata, "Reflection Front")) {
1074	<	if (sd->rb != NULL)
975	<	SDfreeSpectralDF(sd->rb);
976	<	if ((sd->rb = SDnewSpectralDF(1)) == NULL)
1074	>	if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(1)) == NULL)
1075		return SDEmemory;
1076		df = sd->rb;
1077		} else if (!strcasecmp(sdata, "Reflection Back")) {
1078	<	if (sd->rf != NULL)
981	<	SDfreeSpectralDF(sd->rf);
982	<	if ((sd->rf = SDnewSpectralDF(1)) == NULL)
1078	>	if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(1)) == NULL)
1079		return SDEmemory;
1080		df = sd->rf;
1081		} else
1082		return SDEnone;
987	–	/* XXX should also check "ScatteringDataType" for consistency? */
1083		/* get angle bases */
1084		sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis"));
1085		if (!sdata || strcasecmp(sdata, "LBNL/Shirley-Chiu")) {
#	Line 992 | Line 1087 | load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
1087		!sdata ? "Missing" : "Unsupported", sd->name);
1088		return !sdata ? SDEformat : SDEsupport;
1089		}
1090	<	/* allocate BSDF tree */
1091	<	sdt = (SDTre *)malloc(sizeof(SDTre));
1092	<	if (sdt == NULL)
1093	<	return SDEmemory;
1094	<	if (df == sd->rf)
1095	<	sdt->sidef = SD_FREFL;
1096	<	else if (df == sd->rb)
1097	<	sdt->sidef = SD_BREFL;
1098	<	else if (df == sd->tf)
1099	<	sdt->sidef = SD_FXMIT;
1100	<	else /* df == sd->tb */
1101	<	sdt->sidef = SD_BXMIT;
1102	<	sdt->st = NULL;
1103	<	df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
1104	<	df->comp[0].dist = sdt;
1105	<	df->comp[0].func = &SDhandleTre;
1090	>	if (df->comp[0].dist == NULL) { /* need to allocate BSDF tree? */
1091	>	sdt = (SDTre *)malloc(sizeof(SDTre));
1092	>	if (sdt == NULL)
1093	>	return SDEmemory;
1094	>	if (df == sd->rf)
1095	>	sdt->sidef = SD_FREFL;
1096	>	else if (df == sd->rb)
1097	>	sdt->sidef = SD_BREFL;
1098	>	else if (df == sd->tf)
1099	>	sdt->sidef = SD_FXMIT;
1100	>	else /* df == sd->tb */
1101	>	sdt->sidef = SD_BXMIT;
1102	>	sdt->stc[tt_Y] = sdt->stc[tt_u] = sdt->stc[tt_v] = NULL;
1103	>	df->comp[0].dist = sdt;
1104	>	df->comp[0].func = &SDhandleTre;
1105	>	} else {
1106	>	sdt = (SDTre *)df->comp[0].dist;
1107	>	if (sdt->stc[ct] != NULL) {
1108	>	SDfreeTre(sdt->stc[ct]);
1109	>	sdt->stc[ct] = NULL;
1110	>	}
1111	>	}
1112		/* read BSDF data */
1113		sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
1114		if (!sdata || !next_token(&sdata)) {
#	Line 1015 | Line 1116 | load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
1116		sd->name);
1117		return SDEformat;
1118		}
1119	<	sdt->st = load_tree_data(&sdata, ndim);
1120	<	if (sdt->st == NULL)
1119	>	sdt->stc[ct] = load_tree_data(&sdata, ndim);
1120	>	if (sdt->stc[ct] == NULL)
1121		return SDEformat;
1122		if (next_token(&sdata)) {       /* check for unconsumed characters */
1123		sprintf(SDerrorDetail,
#	Line 1025 | Line 1126 | load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
1126		return SDEformat;
1127		}
1128		/* flatten branches where possible */
1129	<	sdt->st = SDsimplifyTre(sdt->st);
1130	<	if (sdt->st == NULL)
1129	>	sdt->stc[ct] = SDsimplifyTre(sdt->stc[ct]);
1130	>	if (sdt->stc[ct] == NULL)
1131		return SDEinternal;
1132	<	return get_extrema(df);         /* compute global quantities */
1132	>	/* compute global quantities for Y */
1133	>	return (ct == tt_Y) ? get_extrema(df) : SDEnone;
1134		}
1135
1136		/* Find minimum value in tree */
#	Line 1068 | Line 1170 | SDsubtractTreVal(SDNode *st, float val)
1170		}
1171		}
1172
1173	<	/* Subtract minimum value from BSDF */
1173	>	/* Subtract minimum Y value from BSDF */
1174		static double
1175	<	subtract_min(SDNode *st)
1175	>	subtract_min_Y(SDNode *st)
1176		{
1177		float   vmin;
1178		/* be sure to skip unused portion */
#	Line 1092 | Line 1194 | subtract_min(SDNode *st)
1194		} else                          /* anisotropic covers entire tree */
1195		vmin = SDgetTreMin(st);
1196
1197	<	if (vmin <= FTINY)
1198	<	return .0;
1197	>	if (vmin <= .01/M_PI)
1198	>	return .0;              /* not worth bothering about */
1199
1200		SDsubtractTreVal(st, vmin);
1201
1202		return M_PI * vmin;             /* return hemispherical value */
1203		}
1204
1205	+	/* Struct used in callback to find RGB extrema */
1206	+	typedef struct {
1207	+	SDNode  **stc;                  /* original Y, u' & v' trees */
1208	+	float   rgb[3];                 /* RGB value */
1209	+	SDNode  *new_stu, *new_stv;     /* replacement u' & v' trees */
1210	+	} SDextRGBs;
1211	+
1212	+	/* Callback to find minimum RGB from Y value plus CIE (u',v') trees */
1213	+	static int
1214	+	get_min_RGB(float yval, const double *cmin, double csiz, void *cptr)
1215	+	{
1216	+	SDextRGBs       *mp = (SDextRGBs *)cptr;
1217	+	double          cmax[SD_MAXDIM];
1218	+	float           rgb[3];
1219	+
1220	+	if (mp->stc[tt_Y]->ndim == 3) {
1221	+	if (cmin[0] + .5*csiz >= .5)
1222	+	return 0;       /* ignore dead half of isotropic */
1223	+	} else
1224	+	cmax[3] = cmin[3] + csiz;
1225	+	cmax[0] = cmin[0] + csiz;
1226	+	cmax[1] = cmin[1] + csiz;
1227	+	cmax[2] = cmin[2] + csiz;
1228	+	/* average RGB color over voxel */
1229	+	SDyuv2rgb(yval, SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1230	+	SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1231	+	/* track smallest components */
1232	+	if (rgb[0] < mp->rgb[0]) mp->rgb[0] = rgb[0];
1233	+	if (rgb[1] < mp->rgb[1]) mp->rgb[1] = rgb[1];
1234	+	if (rgb[2] < mp->rgb[2]) mp->rgb[2] = rgb[2];
1235	+	return 0;
1236	+	}
1237	+
1238	+	/* Callback to build adjusted u' tree */
1239	+	static int
1240	+	adjust_utree(float uprime, const double *cmin, double csiz, void *cptr)
1241	+	{
1242	+	SDextRGBs       *mp = (SDextRGBs *)cptr;
1243	+	double          cmax[SD_MAXDIM];
1244	+	double          yval;
1245	+	float           rgb[3];
1246	+	C_COLOR         clr;
1247	+
1248	+	if (mp->stc[tt_Y]->ndim == 3) {
1249	+	if (cmin[0] + .5*csiz >= .5)
1250	+	return 0;       /* ignore dead half of isotropic */
1251	+	} else
1252	+	cmax[3] = cmin[3] + csiz;
1253	+	cmax[0] = cmin[0] + csiz;
1254	+	cmax[1] = cmin[1] + csiz;
1255	+	cmax[2] = cmin[2] + csiz;
1256	+	/* average RGB color over voxel */
1257	+	SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax), uprime,
1258	+	SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1259	+	/* subtract minimum (& clamp) */
1260	+	if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval;
1261	+	if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval;
1262	+	if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval;
1263	+	c_fromSharpRGB(rgb, &clr);      /* compute new u' for adj. RGB */
1264	+	uprime = 4.*clr.cx/(-2.*clr.cx + 12.*clr.cy + 3.);
1265	+	/* assign in new u' tree */
1266	+	mp->new_stu = SDsetVoxel(mp->new_stu, mp->stc[tt_Y]->ndim,
1267	+	cmin, csiz, uprime);
1268	+	return -(mp->new_stu == NULL);
1269	+	}
1270	+
1271	+	/* Callback to build adjusted v' tree */
1272	+	static int
1273	+	adjust_vtree(float vprime, const double *cmin, double csiz, void *cptr)
1274	+	{
1275	+	SDextRGBs       *mp = (SDextRGBs *)cptr;
1276	+	double          cmax[SD_MAXDIM];
1277	+	double          yval;
1278	+	float           rgb[3];
1279	+	C_COLOR         clr;
1280	+
1281	+	if (mp->stc[tt_Y]->ndim == 3) {
1282	+	if (cmin[0] + .5*csiz >= .5)
1283	+	return 0;       /* ignore dead half of isotropic */
1284	+	} else
1285	+	cmax[3] = cmin[3] + csiz;
1286	+	cmax[0] = cmin[0] + csiz;
1287	+	cmax[1] = cmin[1] + csiz;
1288	+	cmax[2] = cmin[2] + csiz;
1289	+	/* average RGB color over voxel */
1290	+	SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax),
1291	+	SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1292	+	vprime, rgb);
1293	+	/* subtract minimum (& clamp) */
1294	+	if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval;
1295	+	if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval;
1296	+	if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval;
1297	+	c_fromSharpRGB(rgb, &clr);      /* compute new v' for adj. RGB */
1298	+	vprime = 9.*clr.cy/(-2.*clr.cx + 12.*clr.cy + 3.);
1299	+	/* assign in new v' tree */
1300	+	mp->new_stv = SDsetVoxel(mp->new_stv, mp->stc[tt_Y]->ndim,
1301	+	cmin, csiz, vprime);
1302	+	return -(mp->new_stv == NULL);
1303	+	}
1304	+
1305	+	/* Subtract minimum (diffuse) color and return luminance & CIE (x,y) */
1306	+	static double
1307	+	subtract_min_RGB(C_COLOR *cs, SDNode *stc[])
1308	+	{
1309	+	SDextRGBs       my_min;
1310	+	double          ymin;
1311	+
1312	+	my_min.stc = stc;
1313	+	my_min.rgb[0] = my_min.rgb[1] = my_min.rgb[2] = FHUGE;
1314	+	my_min.new_stu = my_min.new_stv = NULL;
1315	+	/* get minimum RGB value */
1316	+	SDtraverseTre(stc[tt_Y], NULL, 0, get_min_RGB, &my_min);
1317	+	/* convert to C_COLOR */
1318	+	ymin =  c_fromSharpRGB(my_min.rgb, cs);
1319	+	if (ymin <= .01/M_PI)           /* not worth bothering about? */
1320	+	return .0;
1321	+	/* adjust u' & v' trees */
1322	+	SDtraverseTre(stc[tt_u], NULL, 0, adjust_utree, &my_min);
1323	+	SDtraverseTre(stc[tt_v], NULL, 0, adjust_vtree, &my_min);
1324	+	SDfreeTre(stc[tt_u]); SDfreeTre(stc[tt_v]);
1325	+	stc[tt_u] = SDsimplifyTre(my_min.new_stu);
1326	+	stc[tt_v] = SDsimplifyTre(my_min.new_stv);
1327	+	/* subtract Y & return hemispherical */
1328	+	SDsubtractTreVal(stc[tt_Y], ymin);
1329	+
1330	+	return M_PI * ymin;
1331	+	}
1332	+
1333		/* Extract and separate diffuse portion of BSDF */
1334		static void
1335		extract_diffuse(SDValue *dv, SDSpectralDF *df)
1336		{
1337		int     n;
1338	+	SDTre   *sdt;
1339
1340		if (df == NULL || df->ncomp <= 0) {
1341		dv->spec = c_dfcolor;
1342		dv->cieY = .0;
1343		return;
1344		}
1345	<	dv->spec = df->comp[0].cspec[0];
1346	<	dv->cieY = subtract_min((*(SDTre *)df->comp[0].dist).st);
1347	<	/* in case of multiple components */
1348	<	for (n = df->ncomp; --n; ) {
1349	<	double  ymin = subtract_min((*(SDTre *)df->comp[n].dist).st);
1350	<	c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
1351	<	dv->cieY += ymin;
1345	>	sdt = (SDTre *)df->comp[0].dist;
1346	>	/* subtract minimum color/grayscale */
1347	>	if (sdt->stc[tt_u] != NULL && sdt->stc[tt_v] != NULL) {
1348	>	int     i = 3*(tt_RGB_coef[1] < .001);
1349	>	while (i--) {           /* initialize on first call */
1350	>	float   rgb[3];
1351	>	rgb[0] = rgb[1] = rgb[2] = .0f; rgb[i] = 1.f;
1352	>	tt_RGB_coef[i] = c_fromSharpRGB(rgb, &tt_RGB_prim[i]);
1353	>	}
1354	>	memcpy(df->comp[0].cspec, tt_RGB_prim, sizeof(tt_RGB_prim));
1355	>	dv->cieY = subtract_min_RGB(&dv->spec, sdt->stc);
1356	>	} else {
1357	>	df->comp[0].cspec[0] = c_dfcolor;
1358	>	dv->cieY = subtract_min_Y(sdt->stc[tt_Y]);
1359		}
1360		df->maxHemi -= dv->cieY;        /* adjust maximum hemispherical */
1361		/* make sure everything is set */
#	Line 1154 | Line 1392 | SDloadTre(SDData *sd, ezxml_t wtl)
1392		/* load BSDF components */
1393		for (wld = ezxml_child(wtl, "WavelengthData");
1394		wld != NULL; wld = wld->next) {
1395	<	if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
1396	<	"Visible"))
1397	<	continue;       /* just visible for now */
1395	>	const char      *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
1396	>	int             ct = -1;
1397	>	if (!strcasecmp(cnm, "Visible"))
1398	>	ct = tt_Y;
1399	>	else if (!strcasecmp(cnm, "CIE-u"))
1400	>	ct = tt_u;
1401	>	else if (!strcasecmp(cnm, "CIE-v"))
1402	>	ct = tt_v;
1403	>	else
1404	>	continue;
1405		for (wdb = ezxml_child(wld, "WavelengthDataBlock");
1406		wdb != NULL; wdb = wdb->next)
1407	<	if ((ec = load_bsdf_data(sd, wdb, rank)) != SDEnone)
1407	>	if ((ec = load_bsdf_data(sd, wdb, ct, rank)) != SDEnone)
1408		return ec;
1409		}
1410		/* separate diffuse components */

Diff Legend

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