--- ray/src/common/bsdf_t.c	2011/04/24 19:39:21	3.6
+++ ray/src/common/bsdf_t.c	2011/04/27 20:03:25	3.7
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf_t.c,v 3.6 2011/04/24 19:39:21 greg Exp $";
+static const char RCSid[] = "$Id: bsdf_t.c,v 3.7 2011/04/27 20:03:25 greg Exp $";
 #endif
 /*
  *  bsdf_t.c
@@ -24,15 +24,18 @@ typedef int	SDtreCallback(float val, const double *cmi
 					double csiz, void *cptr);
 
 					/* reference width maximum (1.0) */
+static const unsigned	iwbits = sizeof(unsigned)*4;
 static const unsigned	iwmax = (1<<(sizeof(unsigned)*4))-1;
+					/* maximum cumulative value */
+static const unsigned	cumlmax = ~0;
 
 /* Struct used for our distribution-building callback */
 typedef struct {
-	int		wmin;		/* minimum square size so far */
-	int		wmax;		/* maximum square size */
 	int		nic;		/* number of input coordinates */
-	int		alen;		/* current array length */
-	int		nall;		/* number of allocated entries */
+	unsigned	alen;		/* current array length */
+	unsigned	nall;		/* number of allocated entries */
+	unsigned	wmin;		/* minimum square size so far */
+	unsigned	wmax;		/* maximum square size */
 	struct outdir_s {
 		unsigned	hent;		/* entering Hilbert index */
 		int		wid;		/* this square size */
@@ -57,12 +60,12 @@ SDnewNode(int nd, int lg)
 	}
 	if (lg < 0) {
 		st = (SDNode *)malloc(sizeof(SDNode) +
-				((1<<nd) - 1)*sizeof(st->u.t[0]));
+				sizeof(st->u.t[0])*((1<<nd) - 1));
 		if (st != NULL)
 			memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
 	} else
 		st = (SDNode *)malloc(sizeof(SDNode) +
-				((1 << nd*lg) - 1)*sizeof(st->u.v[0]));
+				sizeof(st->u.v[0])*((1 << nd*lg) - 1));
 		
 	if (st == NULL) {
 		if (lg < 0)
@@ -103,23 +106,98 @@ SDFreeBTre(void *p)
 	free(sdt);
 }
 
+/* Fill branch's worth of grid values from subtree */
+static void
+fill_grid_branch(float *dptr, const float *sptr, int nd, int shft)
+{
+	unsigned	n = 1 << (shft-1);
+
+	if (!--nd) {			/* end on the line */
+		memcpy(dptr, sptr, sizeof(*dptr)*n);
+		return;
+	}
+	while (n--)			/* recurse on each slice */
+		fill_grid_branch(dptr + (n << shft*nd),
+				sptr + (n << (shft-1)*nd), nd, shft);
+}
+
+/* Get pointer at appropriate offset for the given branch */
+static float *
+grid_branch_start(SDNode *st, int n)
+{
+	unsigned	skipsiz = 1 << st->log2GR;
+	float		*vptr = st->u.v;
+	int		i;
+
+	for (i = st->ndim; i--; skipsiz <<= st->log2GR)
+		if (1<<i & n)
+			vptr += skipsiz >> 1;
+	return vptr;
+}
+
+/* Simplify (consolidate) a tree by flattening uniform depth regions */
+static SDNode *
+SDsimplifyTre(SDNode *st)
+{
+	int		match, n;
+
+	if (st == NULL)			/* check for invalid tree */
+		return NULL;
+	if (st->log2GR >= 0)		/* grid just returns unaltered */
+		return st;
+	match = 1;			/* check if grids below match */
+	for (n = 0; n < 1<<st->ndim; n++) {
+		if ((st->u.t[n] = SDsimplifyTre(st->u.t[n])) == NULL)
+			return NULL;	/* propogate error up call stack */
+		match &= (st->u.t[n]->log2GR == st->u.t[0]->log2GR);
+	}
+	if (match && st->u.t[0]->log2GR >= 0) {
+		SDNode	*stn = SDnewNode(st->ndim, st->u.t[0]->log2GR + 1);
+		if (stn == NULL)	/* out of memory? */
+			return st;
+					/* transfer values to new grid */
+		for (n = 1 << st->ndim; n--; )
+			fill_grid_branch(grid_branch_start(stn, n),
+					st->u.t[n]->u.v, st->ndim, st->log2GR);
+		SDfreeTre(st);		/* free old tree */
+		st = stn;		/* return new one */
+	}
+	return st;
+}
+
+/* Find smallest leaf in tree */
+static double
+SDsmallestLeaf(const SDNode *st)
+{
+	if (st->log2GR < 0) {		/* tree branches */
+		double	lmin = 1.;
+		int	n;
+		for (n = 1<<st->ndim; n--; ) {
+			double	lsiz = SDsmallestLeaf(st->u.t[n]);
+			if (lsiz < lmin)
+				lmin = lsiz;
+		}
+		return .5*lmin;
+	}
+					/* leaf grid width */
+	return 1. / (double)(1 << st->log2GR);
+}
+
 /* Add up N-dimensional hypercube array values over the given box */
 static double
-SDiterSum(const float *va, int nd, int siz, const int *imin, const int *imax)
+SDiterSum(const float *va, int nd, int shft, const int *imin, const int *imax)
 {
+	const unsigned	skipsiz = 1 << nd*shft;
 	double		sum = .0;
-	unsigned	skipsiz = 1;
 	int		i;
 	
-	for (i = nd; --i > 0; )
-		skipsiz *= siz;
 	if (skipsiz == 1)
 		for (i = *imin; i < *imax; i++)
 			sum += va[i];
 	else
 		for (i = *imin; i < *imax; i++)
 			sum += SDiterSum(va + i*skipsiz,
-					nd-1, siz, imin+1, imax+1);
+					nd-1, shft, imin+1, imax+1);
 	return sum;
 }
 
@@ -168,7 +246,7 @@ SDavgTreBox(const SDNode *st, const double *bmin, cons
 	}
 	n = 1;				/* iterate over leaves */
 	for (i = st->ndim; i--; ) {
-		imin[i] = (bmin[i] <= .0) ? 0 
+		imin[i] = (bmin[i] <= 0) ? 0 
 				: (int)((1 << st->log2GR)*bmin[i]);
 		imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR)
 				: (int)((1 << st->log2GR)*bmax[i] + .999999);
@@ -177,8 +255,7 @@ SDavgTreBox(const SDNode *st, const double *bmin, cons
 	if (!n)
 		return .0;
 	
-	return SDiterSum(st->u.v, st->ndim, 1 << st->log2GR, imin, imax) /
-			(double)n;
+	return SDiterSum(st->u.v, st->ndim, st->log2GR, imin, imax) / (double)n;
 }
 
 /* Recursive call for SDtraverseTre() */
@@ -337,9 +414,23 @@ SDqueryTre(const SDTre *sdt, const FVECT outVec, const
 	static const FVECT	zvec = {.0, .0, 1.};
 	FVECT			rOutVec;
 	double			gridPos[4];
-					/* check transmission */
-	if (!sdt->isxmit ^ outVec[2] > 0 ^ inVec[2] > 0)
+
+	switch (sdt->sidef) {		/* whose side are you on? */
+	case SD_UFRONT:
+		if ((outVec[2] < 0) | (inVec[2] < 0))
+			return -1.;
+		break;
+	case SD_UBACK:
+		if ((outVec[2] > 0) | (inVec[2] > 0))
+			return -1.;
+		break;
+	case SD_XMIT:
+		if ((outVec[2] > 0) == (inVec[2] > 0))
+			return -1.;
+		break;
+	default:
 		return -1.;
+	}
 					/* convert vector coordinates */
 	if (sdt->st->ndim == 3) {
 		spinvector(rOutVec, outVec, zvec, -atan2(inVec[1],inVec[0]));
@@ -395,9 +486,8 @@ build_scaffold(float val, const double *cmin, double c
 	bmin[1] = cmin[1]*(double)iwmax + .5;
 	bmax[0] = bmin[0] + wid;
 	bmax[1] = bmin[1] + wid;
-	hilbert_box_vtx(2, sizeof(bitmask_t), sizeof(unsigned)*4,
-							1, bmin, bmax);
-	sp->darr[sp->alen].hent = hilbert_c2i(2, sizeof(unsigned)*4, bmin);
+	hilbert_box_vtx(2, sizeof(bitmask_t), iwbits, 1, bmin, bmax);
+	sp->darr[sp->alen].hent = hilbert_c2i(2, iwbits, bmin);
 	sp->darr[sp->alen].wid = wid;
 	sp->darr[sp->alen].bsdf = val;
 	sp->alen++;			/* on to the next entry */
@@ -416,7 +506,6 @@ sscmp(const void *p1, const void *p2)
 static SDTreCDst *
 make_cdist(const SDTre *sdt, const double *pos)
 {
-	const unsigned	cumlmax = ~0;
 	SDdistScaffold	myScaffold;
 	SDTreCDst	*cd;
 	struct outdir_s	*sp;
@@ -450,14 +539,14 @@ make_cdist(const SDTre *sdt, const double *pos)
 				sizeof(struct outdir_s), &sscmp);
 
 					/* record input range */
-	scale = (double)myScaffold.wmin / iwmax;
+	scale = myScaffold.wmin / (double)iwmax;
 	for (i = myScaffold.nic; i--; ) {
-		cd->clim[i][0] = floor(pos[i]/scale + .5) * scale;
+		cd->clim[i][0] = floor(pos[i]/scale) * scale;
 		cd->clim[i][1] = cd->clim[i][0] + scale;
 	}
 	cd->max_psa = myScaffold.wmax / (double)iwmax;
 	cd->max_psa *= cd->max_psa * M_PI;
-	cd->isxmit = sdt->isxmit;
+	cd->sidef = sdt->sidef;
 	cd->cTotal = 1e-20;		/* compute directional total */
 	sp = myScaffold.darr;
 	for (i = myScaffold.alen; i--; sp++)
@@ -466,6 +555,7 @@ make_cdist(const SDTre *sdt, const double *pos)
 	scale = (double)cumlmax / cd->cTotal;
 	sp = myScaffold.darr;
 	for (i = 0; i < cd->calen; i++, sp++) {
+		cd->carr[i].hndx = sp->hent;
 		cd->carr[i].cuml = scale*cursum + .5;
 		cursum += sp->bsdf * (double)sp->wid * sp->wid;
 	}
@@ -531,9 +621,8 @@ SDqueryTreProjSA(double *psa, const FVECT v1, const RR
 		const SDTre	*sdt = (SDTre *)sdc->dist;
 		double		hcube[SD_MAXDIM];
 		if (SDqueryTre(sdt, v1, v2, hcube) < 0) {
-			if (qflags == SDqueryVal)
-				*psa = M_PI;
-			return SDEnone;
+			strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA");
+			return SDEinternal;
 		}
 		myPSA[0] = hcube[sdt->st->ndim];
 		myPSA[1] = myPSA[0] *= myPSA[0] * M_PI;
@@ -571,15 +660,19 @@ SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
 {
 	const unsigned	nBitsC = 4*sizeof(bitmask_t);
 	const unsigned	nExtraBits = 8*(sizeof(bitmask_t)-sizeof(unsigned));
-	const unsigned	maxval = ~0;
 	const SDTreCDst	*cd = (const SDTreCDst *)cdp;
-	const unsigned	target = randX*maxval;
+	const unsigned	target = randX*cumlmax;
 	bitmask_t	hndx, hcoord[2];
 	double		gpos[3];
 	int		i, iupper, ilower;
 					/* check arguments */
 	if ((ioVec == NULL) | (cd == NULL))
 		return SDEargument;
+	if (ioVec[2] > 0) {
+		if (!(cd->sidef & SD_UFRONT))
+			return SDEargument;
+	} else if (!(cd->sidef & SD_UBACK))
+		return SDEargument;
 					/* binary search to find position */
 	ilower = 0; iupper = cd->calen;
 	while ((i = (iupper + ilower) >> 1) != ilower)
@@ -588,7 +681,7 @@ SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
 		else
 			iupper = i;
 					/* localize random position */
-	randX = (randX*maxval - cd->carr[ilower].cuml) /
+	randX = (randX*cumlmax - cd->carr[ilower].cuml) /
 		    (double)(cd->carr[iupper].cuml - cd->carr[ilower].cuml);
 					/* index in longer Hilbert curve */
 	hndx = (randX*cd->carr[iupper].hndx + (1.-randX)*cd->carr[ilower].hndx)
@@ -599,20 +692,362 @@ SDsampTreCDist(FVECT ioVec, double randX, const SDCDst
 		gpos[i] = ((double)hcoord[i] + rand()*(1./(RAND_MAX+.5))) /
 				(double)((bitmask_t)1 << nBitsC);
 	SDsquare2disk(gpos, gpos[0], gpos[1]);
+					/* compute Z-coordinate */
 	gpos[2] = 1. - gpos[0]*gpos[0] - gpos[1]*gpos[1];
 	if (gpos[2] > 0)		/* paranoia, I hope */
 		gpos[2] = sqrt(gpos[2]);
-	if (ioVec[2] > 0 ^ !cd->isxmit)
+					/* emit from back? */
+	if (ioVec[2] > 0 ^ cd->sidef != SD_XMIT)
 		gpos[2] = -gpos[2];
 	VCOPY(ioVec, gpos);
 	return SDEnone;
 }
 
+/* Advance pointer to the next non-white character in the string (or nul) */
+static int
+next_token(char **spp)
+{
+	while (isspace(**spp))
+		++*spp;
+	return **spp;
+}
+
+/* Count words from this point in string to '}' */
+static int
+count_values(char *cp)
+{
+	int	n = 0;
+
+	while (next_token(&cp) != '}') {
+		if (*cp == '{')
+			return -1;
+		while (*cp && !isspace(*cp))
+			++cp;
+		++n;
+		cp += (next_token(&cp) == ',');
+	}
+	return n;
+}
+
+/* Load an array of real numbers, returning total */
+static int
+load_values(char **spp, float *va, int n)
+{
+	float	*v = va;
+	char	*svnext;
+
+	while (n-- > 0 && (svnext = fskip(*spp)) != NULL) {
+		*v++ = atof(*spp);
+		*spp = svnext;
+		*spp += (next_token(spp) == ',');
+	}
+	return v - va;
+}
+
+/* Load BSDF tree data */
+static SDNode *
+load_tree_data(char **spp, int nd)
+{
+	SDNode	*st;
+	int	n;
+
+	if (next_token(spp) != '{') {
+		strcpy(SDerrorDetail, "Missing '{' in tensor tree");
+		return NULL;
+	}
+	++*spp;				/* in tree, now */
+	if (next_token(spp) == '{') {	/* tree branches */
+		st = SDnewNode(nd, -1);
+		if (st == NULL)
+			return NULL;
+		for (n = 0; n < 1<<nd; n++)
+			if ((st->u.t[n] = load_tree_data(spp, nd)) == NULL) {
+				SDfreeTre(st);
+				return NULL;
+			}
+	} else {			/* else load value grid */
+		int	bsiz;
+		n = count_values(*spp);	/* see how big the grid is */
+		if (n <= 0) {
+			strcpy(SDerrorDetail, "Bad tensor tree data");
+			return NULL;
+		}
+		for (bsiz = 0; bsiz < 8*sizeof(size_t)-1; bsiz += nd)
+			if (1<<bsiz == n)
+				break;
+		if (bsiz >= 8*sizeof(size_t)) {
+			strcpy(SDerrorDetail, "Illegal value count in tensor tree");
+			return NULL;
+		}
+		st = SDnewNode(nd, bsiz/nd);
+		if (st == NULL)
+			return NULL;
+		if (load_values(spp, st->u.v, n) != n) {
+			strcpy(SDerrorDetail, "Real format error in tensor tree");
+			SDfreeTre(st);
+			return NULL;
+		}
+	}
+	if (next_token(spp) != '}') {
+		strcpy(SDerrorDetail, "Missing '}' in tensor tree");
+		SDfreeTre(st);
+		return NULL;
+	}
+	++*spp;				/* walk past close and return */
+	*spp += (next_token(spp) == ',');
+	return st;
+}
+
+/* Compute min. proj. solid angle and max. direct hemispherical scattering */
+static SDError
+get_extrema(SDSpectralDF *df)
+{
+	SDNode	*st = (*(SDTre *)df->comp[0].dist).st;
+	double	stepWidth, dhemi, bmin[4], bmax[4];
+
+	stepWidth = SDsmallestLeaf(st);
+	df->minProjSA = M_PI*stepWidth*stepWidth;
+	if (stepWidth < .03125)
+		stepWidth = .03125;	/* 1/32 resolution good enough */
+	df->maxHemi = .0;
+	if (st->ndim == 3) {		/* isotropic BSDF */
+		bmin[1] = bmin[2] = .0;
+		bmax[1] = bmax[2] = 1.;
+		for (bmin[0] = .0; bmin[0] < .5-FTINY; bmin[0] += stepWidth) {
+			bmax[0] = bmin[0] + stepWidth;
+			dhemi = SDavgTreBox(st, bmin, bmax);
+			if (dhemi > df->maxHemi)
+				df->maxHemi = dhemi;
+		}
+	} else if (st->ndim == 4) {	/* anisotropic BSDF */
+		bmin[2] = bmin[3] = .0;
+		bmax[2] = bmax[3] = 1.;
+		for (bmin[0] = .0; bmin[0] < 1.-FTINY; bmin[0] += stepWidth) {
+			bmax[0] = bmin[0] + stepWidth;
+			for (bmin[1] = .0; bmin[1] < 1.-FTINY; bmin[1] += stepWidth) {
+				bmax[1] = bmin[1] + stepWidth;
+				dhemi = SDavgTreBox(st, bmin, bmax);
+				if (dhemi > df->maxHemi)
+					df->maxHemi = dhemi;
+			}
+		}
+	} else
+		return SDEinternal;
+					/* correct hemispherical value */
+	df->maxHemi *= M_PI;
+	return SDEnone;
+}
+
+/* Load BSDF distribution for this wavelength */
+static SDError
+load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim)
+{
+	SDSpectralDF	*df;
+	SDTre		*sdt;
+	char		*sdata;
+	int		i;
+					/* allocate BSDF component */
+	sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
+	if (!sdata)
+		return SDEnone;
+	/*
+	 * Remember that front and back are reversed from WINDOW 6 orientations
+	 */
+	if (!strcasecmp(sdata, "Transmission")) {
+		if (sd->tf != NULL)
+			SDfreeSpectralDF(sd->tf);
+		if ((sd->tf = SDnewSpectralDF(1)) == NULL)
+			return SDEmemory;
+		df = sd->tf;
+	} else if (!strcasecmp(sdata, "Reflection Front")) {
+		if (sd->rb != NULL)	/* note back-front reversal */
+			SDfreeSpectralDF(sd->rb);
+		if ((sd->rb = SDnewSpectralDF(1)) == NULL)
+			return SDEmemory;
+		df = sd->rb;
+	} else if (!strcasecmp(sdata, "Reflection Back")) {
+		if (sd->rf != NULL)	/* note front-back reversal */
+			SDfreeSpectralDF(sd->rf);
+		if ((sd->rf = SDnewSpectralDF(1)) == NULL)
+			return SDEmemory;
+		df = sd->rf;
+	} else
+		return SDEnone;
+	/* XXX should also check "ScatteringDataType" for consistency? */
+					/* get angle bases */
+	sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis"));
+	if (!sdata || strcasecmp(sdata, "LBNL/Shirley-Chiu")) {
+		sprintf(SDerrorDetail, "%s angle basis for BSDF '%s'",
+				!sdata ? "Missing" : "Unsupported", sd->name);
+		return !sdata ? SDEformat : SDEsupport;
+	}
+					/* allocate BSDF tree */
+	sdt = (SDTre *)malloc(sizeof(SDTre));
+	if (sdt == NULL)
+		return SDEmemory;
+	if (df == sd->rf)
+		sdt->sidef = SD_UFRONT;
+	else if (df == sd->rb)
+		sdt->sidef = SD_UBACK;
+	else
+		sdt->sidef = SD_XMIT;
+	sdt->st = NULL;
+	df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
+	df->comp[0].dist = sdt;
+	df->comp[0].func = &SDhandleTre;
+					/* read BSDF data */
+	sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
+	if (!sdata || !next_token(&sdata)) {
+		sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
+				sd->name);
+		return SDEformat;
+	}
+	sdt->st = load_tree_data(&sdata, ndim);
+	if (sdt->st == NULL)
+		return SDEformat;
+	if (next_token(&sdata)) {	/* check for unconsumed characters */
+		sprintf(SDerrorDetail,
+			"Extra characters at end of ScatteringData in '%s'",
+				sd->name);
+		return SDEformat;
+	}
+					/* flatten branches where possible */
+	sdt->st = SDsimplifyTre(sdt->st);
+	if (sdt->st == NULL)
+		return SDEinternal;
+	return get_extrema(df);		/* compute global quantities */
+}
+
+/* Find minimum value in tree */
+static float
+SDgetTreMin(const SDNode *st)
+{
+	float	vmin = 1./M_PI;
+	int	n;
+
+	if (st->log2GR < 0) {
+		for (n = 1<<st->ndim; n--; ) {
+			float	v = SDgetTreMin(st->u.t[n]);
+			if (v < vmin)
+				vmin = v;
+		}
+	} else {
+		for (n = 1<<(st->ndim*st->log2GR); n--; )
+			if (st->u.v[n] < vmin)
+				vmin = st->u.v[n];
+	}
+	return vmin;
+}
+
+/* Subtract the given value from all tree nodes */
+static void
+SDsubtractTreVal(SDNode *st, float val)
+{
+	int	n;
+
+	if (st->log2GR < 0) {
+		for (n = 1<<st->ndim; n--; )
+			SDsubtractTreVal(st->u.t[n], val);
+	} else {
+		for (n = 1<<(st->ndim*st->log2GR); n--; )
+			st->u.v[n] -= val;
+	}
+}
+
+/* Subtract minimum value from BSDF */
+static double
+subtract_min(SDNode *st)
+{
+	float	vmin;
+					/* be sure to skip unused portion */
+	if ((st->ndim == 3) & (st->log2GR < 0)) {
+		float	v;
+		int	i;
+		vmin = 1./M_PI;
+		for (i = 0; i < 4; i++) {
+			v = SDgetTreMin(st->u.t[i]);
+			if (v < vmin)
+				vmin = v;
+		}
+	} else				/* anisotropic covers entire tree */
+		vmin = SDgetTreMin(st);
+
+	if (vmin <= FTINY)
+		return .0;
+
+	SDsubtractTreMin(st, vmin);
+
+	return M_PI * vmin;		/* return hemispherical value */
+}
+
+/* Extract and separate diffuse portion of BSDF */
+static void
+extract_diffuse(SDValue *dv, SDSpectralDF *df)
+{
+	int	n;
+
+	if (df == NULL || df->ncomp <= 0) {
+		dv->spec = c_dfcolor;
+		dv->cieY = .0;
+		return;
+	}
+	dv->spec = df->comp[0].cspec[0];
+	dv->cieY = subtract_min((*(SDTre *)df->comp[n].dist).st);
+					/* in case of multiple components */
+	for (n = df->ncomp; --n; ) {
+		double	ymin = subtract_min((*(SDTre *)df->comp[n].dist).st);
+		c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
+		dv->cieY += ymin;
+	}
+	df->maxHemi -= dv->cieY;	/* adjust maximum hemispherical */
+					/* make sure everything is set */
+	c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
+}
+
 /* Load a variable-resolution BSDF tree from an open XML file */
 SDError
 SDloadTre(SDData *sd, ezxml_t wtl)
 {
-	return SDEsupport;
+	SDError		ec;
+	ezxml_t		wld, wdb;
+	int		rank;
+	char		*txt;
+					/* basic checks and tensor rank */
+	txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
+			"DataDefinition"), "IncidentDataStructure"));
+	if (txt == NULL || !*txt) {
+		sprintf(SDerrorDetail,
+			"BSDF \"%s\": missing IncidentDataStructure",
+				sd->name);
+		return SDEformat;
+	}
+	if (!strcasecmp(txt, "TensorTree3"))
+		rank = 3;
+	else if (!strcasecmp(txt, "TensorTree4"))
+		rank = 4;
+	else {
+		sprintf(SDerrorDetail,
+			"BSDF \"%s\": unsupported IncidentDataStructure",
+				sd->name);
+		return SDEsupport;
+	}
+					/* load BSDF components */
+	for (wld = ezxml_child(wtl, "WavelengthData");
+				wld != NULL; wld = wld->next) {
+		if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
+				"Visible"))
+			continue;	/* just visible for now */
+		for (wdb = ezxml_child(wld, "WavelengthDataBlock");
+					wdb != NULL; wdb = wdb->next)
+			if ((ec = load_bsdf_data(sd, wdb, rank)) != SDEnone)
+				return ec;
+	}
+					/* separate diffuse components */
+	extract_diffuse(&sd->rLambFront, sd->rf);
+	extract_diffuse(&sd->rLambBack, sd->rb);
+	extract_diffuse(&sd->tLamb, sd->tf);
+					/* return success */
+	return SDEnone;
 }
 
 /* Variable resolution BSDF methods */