ViewVC Help

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

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

Comparing ray/src/common/bsdf.c (file contents):
Revision 2.18 by greg, Fri Apr 8 23:23:28 2011 UTC vs.
Revision 2.37 by greg, Sun Mar 4 20:11:10 2012 UTC

#	Line 10 | Line 10 | static const char RCSid[] = "$Id$";
10		*
11		*/
12
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"
#	Line 45 | Line 48 | int                    SDretainSet = SDretainNone;
48		SDError
49		SDreportEnglish(SDError ec, FILE *fp)
50		{
48	–	if (fp == NULL)
49	–	return ec;
51		if (!ec)
52		return SDEnone;
53	+	if ((ec < SDEnone) | (ec > SDEunknown)) {
54	+	SDerrorDetail[0] = '\0';
55	+	ec = SDEunknown;
56	+	}
57	+	if (fp == NULL)
58	+	return ec;
59		fputs(SDerrorEnglish[ec], fp);
60		if (SDerrorDetail[0]) {
61		fputs(": ", fp);
#	Line 95 | Line 102 | SDloadGeometry(SDData *sd, ezxml_t wdb)
102		if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
103		sd->dim[2] = atof(ezxml_txt(geom)) *
104		to_meters(ezxml_attr(geom, "unit"));
105	<	if ((sd->dim[0] < .0) | (sd->dim[1] < .0) | (sd->dim[2] < .0)) {
105	>	if ((sd->dim[0] < 0) | (sd->dim[1] < 0) | (sd->dim[2] < 0)) {
106		sprintf(SDerrorDetail, "Negative size in \"%s\"", sd->name);
107		return SDEdata;
108		}
109		if ((geom = ezxml_child(wdb, "Geometry")) == NULL ||
110		(mgfstr = ezxml_txt(geom)) == NULL)
111		return SDEnone;
112	+	while (isspace(*mgfstr))
113	+	++mgfstr;
114	+	if (!*mgfstr)
115	+	return SDEnone;
116		if ((fmt = ezxml_attr(geom, "format")) != NULL &&
117		strcasecmp(fmt, "MGF")) {
118		sprintf(SDerrorDetail,
#	Line 160 | Line 171 | SDloadFile(SDData *sd, const char *fname)
171		}
172		/* load geometry if present */
173		lastErr = SDloadGeometry(sd, ezxml_child(wtl, "Material"));
174	<	if (lastErr)
174	>	if (lastErr) {
175	>	ezxml_free(fl);
176		return lastErr;
177	+	}
178		/* try loading variable resolution data */
179		lastErr = SDloadTre(sd, wtl);
180		/* check our result */
181	<	switch (lastErr) {
169	<	case SDEformat:
170	<	case SDEdata:
171	<	case SDEsupport:        /* possibly we just tried the wrong format */
181	>	if (lastErr == SDEsupport)      /* try matrix BSDF if not tree data */
182		lastErr = SDloadMtx(sd, wtl);
183	<	break;
174	<	default:                /* variable res. OK else serious error */
175	<	break;
176	<	}
183	>
184		/* done with XML file */
185		ezxml_free(fl);
186
#	Line 219 | Line 226 | SDnewSpectralDF(int nc)
226		return df;
227		}
228
229	+	/* Add component(s) to spectral distribution function */
230	+	SDSpectralDF *
231	+	SDaddComponent(SDSpectralDF *odf, int nadd)
232	+	{
233	+	SDSpectralDF    *df;
234	+
235	+	if (odf == NULL)
236	+	return SDnewSpectralDF(nadd);
237	+	if (nadd <= 0)
238	+	return odf;
239	+	df = (SDSpectralDF *)realloc(odf, sizeof(SDSpectralDF) +
240	+	(odf->ncomp+nadd-1)*sizeof(SDComponent));
241	+	if (df == NULL) {
242	+	sprintf(SDerrorDetail,
243	+	"Cannot add %d component(s) to spectral DF", nadd);
244	+	SDfreeSpectralDF(odf);
245	+	return NULL;
246	+	}
247	+	memset(df->comp+df->ncomp, 0, nadd*sizeof(SDComponent));
248	+	df->ncomp += nadd;
249	+	return df;
250	+	}
251	+
252		/* Free cached cumulative distributions for BSDF component */
253		void
254		SDfreeCumulativeCache(SDSpectralDF *df)
#	Line 245 | Line 275 | SDfreeSpectralDF(SDSpectralDF *df)
275		return;
276		SDfreeCumulativeCache(df);
277		for (n = df->ncomp; n-- > 0; )
278	<	(*df->comp[n].func->freeSC)(df->comp[n].dist);
278	>	if (df->comp[n].dist != NULL)
279	>	(*df->comp[n].func->freeSC)(df->comp[n].dist);
280		free(df);
281		}
282
#	Line 269 | Line 300 | SDclipName(char *res, const char *fname)
300
301		/* Initialize an unused BSDF struct (simply clears to zeroes) */
302		void
303	<	SDclearBSDF(SDData *sd)
303	>	SDclearBSDF(SDData *sd, const char *fname)
304		{
305	<	if (sd != NULL)
306	<	memset(sd, 0, sizeof(SDData));
305	>	if (sd == NULL)
306	>	return;
307	>	memset(sd, 0, sizeof(SDData));
308	>	if (fname == NULL)
309	>	return;
310	>	SDclipName(sd->name, fname);
311		}
312
313		/* Free data associated with BSDF struct */
#	Line 330 | Line 365 | SDgetCache(const char *bname)
365		sdl->next = SDcacheList;
366		SDcacheList = sdl;
367
368	<	sdl->refcnt++;
368	>	sdl->refcnt = 1;
369		return &sdl->bsdf;
370		}
371
#	Line 374 | Line 409 | SDfreeCache(const SDData *sd)
409		for (sdl = SDcacheList; sdl != NULL; sdl = (sdLast=sdl)->next)
410		if (&sdl->bsdf == sd)
411		break;
412	<	if (sdl == NULL || --sdl->refcnt)
412	>	if (sdl == NULL || (sdl->refcnt -= (sdl->refcnt > 0)))
413		return;                 /* missing or still in use */
414		/* keep unreferenced data? */
415		if (SDisLoaded(sd) && SDretainSet) {
#	Line 397 | Line 432 | SDfreeCache(const SDData *sd)
432
433		/* Sample an individual BSDF component */
434		SDError
435	<	SDsampComponent(SDValue *sv, FVECT outVec, const FVECT inVec,
401	<	double randX, SDComponent *sdc)
435	>	SDsampComponent(SDValue *sv, FVECT ioVec, double randX, SDComponent *sdc)
436		{
437		float           coef[SDmaxCh];
438		SDError         ec;
439	+	FVECT           inVec;
440		const SDCDst    *cd;
441		double          d;
442		int             n;
443		/* check arguments */
444	<	if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL))
444	>	if ((sv == NULL) | (ioVec == NULL) | (sdc == NULL))
445		return SDEargument;
446		/* get cumulative distribution */
447	+	VCOPY(inVec, ioVec);
448		cd = (*sdc->func->getCDist)(inVec, sdc);
449		if (cd == NULL)
450		return SDEmemory;
451	<	if (cd->cTotal <= 1e-7) {       /* anything to sample? */
451	>	if (cd->cTotal <= 1e-6) {       /* anything to sample? */
452		sv->spec = c_dfcolor;
453		sv->cieY = .0;
454	<	memset(outVec, 0, 3*sizeof(double));
454	>	memset(ioVec, 0, 3*sizeof(double));
455		return SDEnone;
456		}
457		sv->cieY = cd->cTotal;
458		/* compute sample direction */
459	<	ec = (*sdc->func->sampCDist)(outVec, randX, cd);
459	>	ec = (*sdc->func->sampCDist)(ioVec, randX, cd);
460		if (ec)
461		return ec;
462		/* get BSDF color */
463	<	n = (*sdc->func->getBSDFs)(coef, outVec, inVec, sdc->dist);
463	>	n = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
464		if (n <= 0) {
465		strcpy(SDerrorDetail, "BSDF sample value error");
466		return SDEinternal;
#	Line 451 | Line 487 | SDmultiSamp(double t[], int n, double randX)
487		bitmask_t       ndx, coord[MS_MAXDIM];
488
489		while (n > MS_MAXDIM)           /* punt for higher dimensions */
490	<	t[--n] = rand()*(1./RAND_MAX);
490	>	t[--n] = rand()*(1./(RAND_MAX+.5));
491		nBits = (8*sizeof(bitmask_t) - 1) / n;
492		ndx = randX * (double)((bitmask_t)1 << (nBits*n));
493		/* get coordinate on Hilbert curve */
#	Line 459 | Line 495 | SDmultiSamp(double t[], int n, double randX)
495		/* convert back to [0,1) range */
496		scale = 1. / (double)((bitmask_t)1 << nBits);
497		while (n--)
498	<	t[n] = scale * ((double)coord[n] + rand()*(1./RAND_MAX));
498	>	t[n] = scale * ((double)coord[n] + rand()*(1./(RAND_MAX+.5)));
499		}
500
501		#undef MS_MAXDIM
#	Line 472 | Line 508 | SDdiffuseSamp(FVECT outVec, int outFront, double randX
508		SDmultiSamp(outVec, 2, randX);
509		SDsquare2disk(outVec, outVec[0], outVec[1]);
510		outVec[2] = 1. - outVec[0]*outVec[0] - outVec[1]*outVec[1];
511	<	if (outVec[2] > .0)             /* a bit of paranoia */
511	>	if (outVec[2] > 0)              /* a bit of paranoia */
512		outVec[2] = sqrt(outVec[2]);
513		if (!outFront)                  /* going out back? */
514		outVec[2] = -outVec[2];
#	Line 480 | Line 516 | SDdiffuseSamp(FVECT outVec, int outFront, double randX
516
517		/* Query projected solid angle coverage for non-diffuse BSDF direction */
518		SDError
519	<	SDsizeBSDF(double *projSA, const FVECT vec, int qflags, const SDData *sd)
519	>	SDsizeBSDF(double *projSA, const FVECT v1, const RREAL *v2,
520	>	int qflags, const SDData *sd)
521		{
522	<	SDSpectralDF    *rdf;
522	>	SDSpectralDF    *rdf, *tdf;
523		SDError         ec;
524		int             i;
525		/* check arguments */
526	<	if ((projSA == NULL) | (vec == NULL) | (sd == NULL))
526	>	if ((projSA == NULL) | (v1 == NULL) | (sd == NULL))
527		return SDEargument;
528		/* initialize extrema */
529		switch (qflags) {
#	Line 502 | Line 539 | SDsizeBSDF(double *projSA, const FVECT vec, int qflags
539		case 0:
540		return SDEargument;
541		}
542	<	if (vec[2] > .0)                /* front surface query? */
542	>	if (v1[2] > 0)                  /* front surface query? */
543		rdf = sd->rf;
544		else
545		rdf = sd->rb;
546	+	tdf = sd->tf;
547	+	if (v2 != NULL)                 /* bidirectional? */
548	+	if (v1[2] > 0 ^ v2[2] > 0)
549	+	rdf = NULL;
550	+	else
551	+	tdf = NULL;
552		ec = SDEdata;                   /* run through components */
553		for (i = (rdf==NULL) ? 0 : rdf->ncomp; i--; ) {
554	<	ec = (*rdf->comp[i].func->queryProjSA)(projSA, vec, qflags,
555	<	rdf->comp[i].dist);
554	>	ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2,
555	>	qflags, &rdf->comp[i]);
556		if (ec)
557		return ec;
558		}
559	<	for (i = (sd->tf==NULL) ? 0 : sd->tf->ncomp; i--; ) {
560	<	ec = (*sd->tf->comp[i].func->queryProjSA)(projSA, vec, qflags,
561	<	sd->tf->comp[i].dist);
559	>	for (i = (tdf==NULL) ? 0 : tdf->ncomp; i--; ) {
560	>	ec = (*tdf->comp[i].func->queryProjSA)(projSA, v1, v2,
561	>	qflags, &tdf->comp[i]);
562		if (ec)
563		return ec;
564		}
#	Line 539 | Line 582 | SDevalBSDF(SDValue *sv, const FVECT outVec, const FVEC
582		if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL))
583		return SDEargument;
584		/* whose side are we on? */
585	<	inFront = (inVec[2] > .0);
586	<	outFront = (outVec[2] > .0);
585	>	inFront = (inVec[2] > 0);
586	>	outFront = (outVec[2] > 0);
587		/* start with diffuse portion */
588		if (inFront & outFront) {
589		*sv = sd->rLambFront;
#	Line 557 | Line 600 | SDevalBSDF(SDValue *sv, const FVECT outVec, const FVEC
600		i = (sdf != NULL) ? sdf->ncomp : 0;
601		while (i-- > 0) {
602		nch = (*sdf->comp[i].func->getBSDFs)(coef, outVec, inVec,
603	<	sdf->comp[i].dist);
603	>	&sdf->comp[i]);
604		while (nch-- > 0) {
605		c_cmix(&sv->spec, sv->cieY, &sv->spec,
606		coef[nch], &sdf->comp[i].cspec[nch]);
#	Line 581 | Line 624 | SDdirectHemi(const FVECT inVec, int sflags, const SDDa
624		if ((inVec == NULL) | (sd == NULL))
625		return .0;
626		/* gather diffuse components */
627	<	if (inVec[2] > .0) {
627	>	if (inVec[2] > 0) {
628		hsum = sd->rLambFront.cieY;
629		rdf = sd->rf;
630		} else /* !inFront */ {
#	Line 612 | Line 655 | SDdirectHemi(const FVECT inVec, int sflags, const SDDa
655
656		/* Sample BSDF direction based on the given random variable */
657		SDError
658	<	SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVec,
616	<	double randX, int sflags, const SDData *sd)
658	>	SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int sflags, const SDData *sd)
659		{
660		SDError         ec;
661	+	FVECT           inVec;
662		int             inFront;
663		SDSpectralDF    *rdf;
664		double          rdiff;
#	Line 624 | Line 667 | SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
667		SDComponent     *sdc;
668		const SDCDst    **cdarr = NULL;
669		/* check arguments */
670	<	if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL) |
671	<	(randX < .0) | (randX >= 1.))
670	>	if ((sv == NULL) | (ioVec == NULL) | (sd == NULL) |
671	>	(randX < 0) | (randX >= 1.))
672		return SDEargument;
673		/* whose side are we on? */
674	<	inFront = (inVec[2] > .0);
674	>	VCOPY(inVec, ioVec);
675	>	inFront = (inVec[2] > 0);
676		/* remember diffuse portions */
677		if (inFront) {
678		*sv = sd->rLambFront;
#	Line 666 | Line 710 | SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
710		}
711		sv->cieY += cdarr[i]->cTotal;
712		}
713	<	if (sv->cieY <= 1e-7) {         /* anything to sample? */
713	>	if (sv->cieY <= 1e-6) {         /* anything to sample? */
714		sv->cieY = .0;
715	<	memset(outVec, 0, 3*sizeof(double));
715	>	memset(ioVec, 0, 3*sizeof(double));
716		return SDEnone;
717		}
718		/* scale random variable */
719		randX *= sv->cieY;
720		/* diffuse reflection? */
721		if (randX < rdiff) {
722	<	SDdiffuseSamp(outVec, inFront, randX/rdiff);
722	>	SDdiffuseSamp(ioVec, inFront, randX/rdiff);
723		goto done;
724		}
725		randX -= rdiff;
#	Line 683 | Line 727 | SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
727		if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) {
728		if (randX < sd->tLamb.cieY) {
729		sv->spec = sd->tLamb.spec;
730	<	SDdiffuseSamp(outVec, !inFront, randX/sd->tLamb.cieY);
730	>	SDdiffuseSamp(ioVec, !inFront, randX/sd->tLamb.cieY);
731		goto done;
732		}
733		randX -= sd->tLamb.cieY;
#	Line 695 | Line 739 | SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe
739		return SDEinternal;
740		/* compute sample direction */
741		sdc = (i < nr) ? &rdf->comp[i] : &sd->tf->comp[i-nr];
742	<	ec = (*sdc->func->sampCDist)(outVec, randX/cdarr[i]->cTotal, cdarr[i]);
742	>	ec = (*sdc->func->sampCDist)(ioVec, randX/cdarr[i]->cTotal, cdarr[i]);
743		if (ec)
744		return ec;
745		/* compute color */
746	<	j = (*sdc->func->getBSDFs)(coef, outVec, inVec, sdc->dist);
746	>	j = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
747		if (j <= 0) {
748		sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error",
749		sd->name);
#	Line 726 | Line 770 | SDcompXform(RREAL vMtx[3][3], const FVECT sNrm, const
770		if ((vMtx == NULL) | (sNrm == NULL) | (uVec == NULL))
771		return SDEargument;
772		VCOPY(vMtx[2], sNrm);
773	<	if (normalize(vMtx[2]) == .0)
773	>	if (normalize(vMtx[2]) == 0)
774		return SDEargument;
775		fcross(vMtx[0], uVec, vMtx[2]);
776	<	if (normalize(vMtx[0]) == .0)
776	>	if (normalize(vMtx[0]) == 0)
777		return SDEargument;
778		fcross(vMtx[1], vMtx[2], vMtx[0]);
779		return SDEnone;
#	Line 749 | Line 793 | SDinvXform(RREAL iMtx[3][3], RREAL vMtx[3][3])
793		mTmp[0][1] = vMtx[2][1]*vMtx[0][2] - vMtx[2][2]*vMtx[0][1];
794		mTmp[0][2] = vMtx[1][2]*vMtx[0][1] - vMtx[1][1]*vMtx[0][2];
795		d = vMtx[0][0]*mTmp[0][0] + vMtx[1][0]*mTmp[0][1] + vMtx[2][0]*mTmp[0][2];
796	<	if (d == .0) {
796	>	if (d == 0) {
797		strcpy(SDerrorDetail, "Zero determinant in matrix inversion");
798		return SDEargument;
799		}
#	Line 776 | Line 820 | SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT i
820		if (vMtx == NULL) {             /* assume they just want to normalize */
821		if (resVec != inpVec)
822		VCOPY(resVec, inpVec);
823	<	return (normalize(resVec) > .0) ? SDEnone : SDEargument;
823	>	return (normalize(resVec) > 0) ? SDEnone : SDEargument;
824		}
825		vTmp[0] = DOT(vMtx[0], inpVec);
826		vTmp[1] = DOT(vMtx[1], inpVec);
827		vTmp[2] = DOT(vMtx[2], inpVec);
828	<	if (normalize(vTmp) == .0)
828	>	if (normalize(vTmp) == 0)
829		return SDEargument;
830		VCOPY(resVec, vTmp);
831		return SDEnone;
#	Line 796 | Line 840 | SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT i
840
841		#include "standard.h"
842		#include "paths.h"
799	–	#include <ctype.h>
843
844		#define MAXLATS         46              /* maximum number of latitudes */
845
#	Line 853 | Line 896 | static int     nabases = 3;    /* current number of defined b
896		static int
897		fequal(double a, double b)
898		{
899	<	if (b != .0)
899	>	if (b != 0)
900		a = a/b - 1.;
901		return((a <= 1e-6) & (a >= -1e-6));
902		}
#	Line 916 | Line 959 | ab_getndx(             /* get index corresponding to the given ve
959		{
960		ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
961		int     li, ndx;
962	<	double  pol, azi, d;
962	>	double  pol, azi;
963
964		if ((v[2] < -1.0) | (v[2] > 1.0))
965		return(-1);
#	Line 1156 | Line 1199 | check_bsdf_data(       /* check that BSDF data is sane */
1199		)
1200		{
1201		double          *omega_iarr, *omega_oarr;
1202	<	double          dom, contrib, hemi_total, full_total;
1202	>	double          dom, hemi_total, full_total;
1203		int             nneg;
1204		FVECT           v;
1205		int             i, o;
#	Line 1171 | Line 1214 | check_bsdf_data(       /* check that BSDF data is sane */
1214		hemi_total = .0;
1215		for (i = dp->ninc; i--; ) {
1216		dom = getBSDF_incohm(dp,i);
1217	<	if (dom <= .0) {
1217	>	if (dom <= 0) {
1218		error(WARNING, "zero/negative incoming solid angle");
1219		continue;
1220		}
#	Line 1194 | Line 1237 | check_bsdf_data(       /* check that BSDF data is sane */
1237		hemi_total = .0;
1238		for (o = dp->nout; o--; ) {
1239		dom = getBSDF_outohm(dp,o);
1240	<	if (dom <= .0) {
1240	>	if (dom <= 0) {
1241		error(WARNING, "zero/negative outgoing solid angle");
1242		continue;
1243		}
#	Line 1218 | Line 1261 | check_bsdf_data(       /* check that BSDF data is sane */
1261		hemi_total = .0;
1262		for (o = dp->nout; o--; ) {
1263		double  f = BSDF_value(dp,i,o);
1264	<	if (f >= .0)
1264	>	if (f >= 0)
1265		hemi_total += f*omega_oarr[o];
1266		else {
1267		nneg += (f < -FTINY);
#	Line 1306 | Line 1349 | load_BSDF(             /* load BSDF data from file */
1349		error(WARNING, errmsg);
1350		ezxml_free(fl);
1351		return(NULL);
1352	<	}
1353	<	load_angle_basis(ezxml_child(ezxml_child(wtl,
1354	<	"DataDefinition"), "AngleBasis"));
1352	>	}
1353	>	for (wld = ezxml_child(ezxml_child(wtl,
1354	>	"DataDefinition"), "AngleBasis");
1355	>	wld != NULL; wld = wld->next)
1356	>	load_angle_basis(wld);
1357		dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
1358		load_geometry(dp, ezxml_child(wtl, "Material"));
1359		for (wld = ezxml_child(wtl, "WavelengthData");

Diff Legend

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