--- ray/src/common/bsdf_m.c	2015/04/05 01:32:01	3.30
+++ ray/src/common/bsdf_m.c	2022/01/25 01:34:20	3.45
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf_m.c,v 3.30 2015/04/05 01:32:01 greg Exp $";
+static const char RCSid[] = "$Id: bsdf_m.c,v 3.45 2022/01/25 01:34:20 greg Exp $";
 #endif
 /*
  *  bsdf_m.c
@@ -132,6 +132,9 @@ SDfreeMatrix(void *ptr)
 	free(ptr);
 }
 
+/* compute square of real value */
+static double sq(double x) { return x*x; }
+
 /* Get vector for this angle basis index (front exiting) */
 int
 fo_getvec(FVECT v, double ndxr, void *p)
@@ -141,17 +144,17 @@ fo_getvec(FVECT v, double ndxr, void *p)
 	double		randX = ndxr - ndx;
 	double		rx[2];
 	int		li;
-	double		pol, azi, d;
+	double		azi, d;
 	
 	if ((ndxr < 0) | (ndx >= ab->nangles))
 		return RC_FAIL;
 	for (li = 0; ndx >= ab->lat[li].nphis; li++)
 		ndx -= ab->lat[li].nphis;
 	SDmultiSamp(rx, 2, randX);
-	pol = M_PI/180.*( (1.-rx[0])*ab->lat[li].tmin +
-				rx[0]*ab->lat[li+1].tmin );
+	d = (1. - rx[0])*sq(cos(M_PI/180.*ab->lat[li].tmin)) +
+		rx[0]*sq(cos(M_PI/180.*ab->lat[li+1].tmin));
+	v[2] = d = sqrt(d);	/* cos(pol) */
 	azi = 2.*M_PI*(ndx + rx[1] - .5)/ab->lat[li].nphis;
-	v[2] = d = cos(pol);
 	d = sqrt(1. - d*d);	/* sin(pol) */
 	v[0] = cos(azi)*d;
 	v[1] = sin(azi)*d;
@@ -168,7 +171,7 @@ fo_getndx(const FVECT v, void *p)
 
 	if (v == NULL)
 		return -1;
-	if ((v[2] < 0) | (v[2] > 1.))
+	if ((v[2] < 0) | (v[2] > 1.00001))
 		return -1;
 	pol = 180.0/M_PI*Acos(v[2]);
 	azi = 180.0/M_PI*atan2(v[1], v[0]);
@@ -184,13 +187,11 @@ fo_getndx(const FVECT v, void *p)
 	return ndx;
 }
 
-/* compute square of real value */
-static double sq(double x) { return x*x; }
-
 /* Get projected solid angle for this angle basis index (universal) */
 double
 io_getohm(int ndx, void *p)
 {
+	static void	*last_p = NULL;
 	static int	last_li = -1;
 	static double	last_ohm;
 	ANGLE_BASIS	*ab = (ANGLE_BASIS *)p;
@@ -201,8 +202,9 @@ io_getohm(int ndx, void *p)
 		return -1.;
 	for (li = 0; ndx >= ab->lat[li].nphis; li++)
 		ndx -= ab->lat[li].nphis;
-	if (li == last_li)			/* cached latitude? */
+	if ((p == last_p) & (li == last_li))		/* cached latitude? */
 		return last_ohm;
+	last_p = p;
 	last_li = li;
 	theta = M_PI/180. * ab->lat[li].tmin;
 	theta1 = M_PI/180. * ab->lat[li+1].tmin;
@@ -293,12 +295,17 @@ int
 mBSDF_color(float coef[], const SDMat *dp, int i, int o)
 {
 	C_COLOR	cxy;
+	double	d;
 
-	coef[0] = mBSDF_value(dp, i, o);
+	coef[0] = mBSDF_value(dp, o, i);
+				/* position-specific perturbation */
+	d = 2*dp->ninc/(i + .22545) + 4*dp->nout/(o + .70281);
+	d -= (int)d;
+	coef[0] *= 1. + 6e-4*(d - .5);
 	if (dp->chroma == NULL)
 		return 1;	/* grayscale */
 
-	c_decodeChroma(&cxy, dp->chroma[o*dp->ninc + i]);
+	c_decodeChroma(&cxy, mBSDF_chroma(dp,o,i));
 	c_toSharpRGB(&cxy, coef[0], coef);
 	coef[0] *= mtx_RGB_coef[0];
 	coef[1] *= mtx_RGB_coef[1];
@@ -381,7 +388,7 @@ get_extrema(SDSpectralDF *df)
 	for (i = dp->ninc; i--; ) {
 		double	hemi = .0;
 		for (o = dp->nout; o--; )
-			hemi += ohma[o] * mBSDF_value(dp, i, o);
+			hemi += ohma[o] * mBSDF_value(dp, o, i);
 		if (hemi > df->maxHemi)
 			df->maxHemi = hemi;
 	}
@@ -517,7 +524,7 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int ro
 		if (rowinc) {
 			int	r = i/dp->nout;
 			int	c = i - r*dp->nout;
-			mBSDF_value(dp,r,c) = val;
+			mBSDF_value(dp,c,r) = val;
 		} else
 			dp->bsdf[i] = val;
 		sdata = sdnext;
@@ -591,7 +598,7 @@ static double
 subtract_min(C_COLOR *cs, SDMat *sm)
 {
 	const int	ncomp = 1 + 2*(sm->chroma != NULL);
-	float		min_coef[3], coef[3];
+	float		min_coef[3], ymin, coef[3];
 	int		i, o, c;
 	
 	min_coef[0] = min_coef[1] = min_coef[2] = FHUGE;
@@ -602,16 +609,16 @@ subtract_min(C_COLOR *cs, SDMat *sm)
 				if (coef[c] < min_coef[c])
 					min_coef[c] = coef[c];
 		}
+	ymin = 0;
 	for (c = ncomp; c--; )
-		if (min_coef[c] > FTINY)
-			break;
-	if (c < 0)
+		ymin += min_coef[c];
+	if (ymin <= .01/M_PI)		/* not worth bothering about? */
 		return .0;
 	if (ncomp == 1) {		/* subtract grayscale minimum */
 		for (i = sm->ninc*sm->nout; i--; )
-			sm->bsdf[i] -= min_coef[0];
+			sm->bsdf[i] -= ymin;
 		*cs = c_dfcolor;
-		return min_coef[0]*M_PI;
+		return M_PI*ymin;
 	}
 					/* else subtract colored minimum */
 	for (i = 0; i < sm->ninc; i++)
@@ -622,15 +629,15 @@ subtract_min(C_COLOR *cs, SDMat *sm)
 				coef[c] = (coef[c] - min_coef[c]) /
 						mtx_RGB_coef[c];
 			if (c_fromSharpRGB(coef, &cxy) > 1e-5)
-				sm->chroma[o*sm->ninc + i] = c_encodeChroma(&cxy);
-			mBSDF_value(sm,i,o) -= min_coef[0]+min_coef[1]+min_coef[2];
+				mBSDF_chroma(sm,o,i) = c_encodeChroma(&cxy);
+			mBSDF_value(sm,o,i) -= ymin;
 		}
 					/* return colored minimum */
 	for (i = 3; i--; )
 		coef[i] = min_coef[i]/mtx_RGB_coef[i];
 	c_fromSharpRGB(coef, cs);
 
-	return (min_coef[0]+min_coef[1]+min_coef[2])*M_PI;
+	return M_PI*ymin;
 }
 
 /* Extract and separate diffuse portion of BSDF & convert color */
@@ -647,8 +654,8 @@ extract_diffuse(SDValue *dv, SDSpectralDF *df)
 					/* subtract minimum value */
 	dv->cieY = subtract_min(&dv->spec, (SDMat *)df->comp[0].dist);
 	df->maxHemi -= dv->cieY;	/* adjust maximum hemispherical */
-					/* make sure everything is set */
-	c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
+				
+	c_ccvt(&dv->spec, C_CSXY);	/* make sure (x,y) is set */
 	return df;
 }
 
@@ -707,18 +714,21 @@ SDloadMtx(SDData *sd, ezxml_t wtl)
 					/* separate diffuse components */
 	sd->rf = extract_diffuse(&sd->rLambFront, sd->rf);
 	sd->rb = extract_diffuse(&sd->rLambBack, sd->rb);
-	if (sd->tf != NULL)
-		sd->tf = extract_diffuse(&sd->tLamb, sd->tf);
-	if (sd->tb != NULL)
-		sd->tb = extract_diffuse(&sd->tLamb, sd->tb);
+	sd->tf = extract_diffuse(&sd->tLambFront, sd->tf);
+	if (sd->tb != NULL) {
+		sd->tb = extract_diffuse(&sd->tLambBack, sd->tb);
+		if (sd->tf == NULL)
+			sd->tLambFront = sd->tLambBack;
+	} else if (sd->tf != NULL)
+		sd->tLambBack = sd->tLambFront;
 					/* return success */
 	return SDEnone;
 }
 
 /* Get Matrix BSDF value */
 static int
-SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
-				const FVECT inVec, SDComponent *sdc)
+SDgetMtxBSDF(float coef[SDmaxCh], const FVECT inVec,
+				const FVECT outVec, SDComponent *sdc)
 {
 	const SDMat	*dp;
 	int		i_ndx, o_ndx;
@@ -803,10 +813,10 @@ make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp
 	cmtab[0] = .0;
 	for (o = 0; o < cd->calen; o++) {
 		if (rev)
-			cmtab[o+1] = mBSDF_value(dp, o, cd->indx) *
+			cmtab[o+1] = mBSDF_value(dp, cd->indx, o) *
 					(*dp->ib_ohm)(o, dp->ib_priv);
 		else
-			cmtab[o+1] = mBSDF_value(dp, cd->indx, o) *
+			cmtab[o+1] = mBSDF_value(dp, o, cd->indx) *
 					(*dp->ob_ohm)(o, dp->ob_priv);
 		cmtab[o+1] += cmtab[o];
 	}
@@ -849,6 +859,7 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 		reverse = 1;
 	}
 	cdlast = NULL;			/* check for it in cache list */
+	/* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
 	for (cd = (SDMatCDst *)sdc->cdList; cd != NULL;
 					cdlast = cd, cd = cd->next)
 		if (cd->indx == myCD.indx && (cd->calen == myCD.calen) &
@@ -872,6 +883,7 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 		cd->next = (SDMatCDst *)sdc->cdList;
 		sdc->cdList = (SDCDst *)cd;
 	}
+	/* END MUTEX LOCK */
 	return (SDCDst *)cd;		/* ready to go */
 }
 
@@ -904,7 +916,7 @@ SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst
 }
 
 /* Fixed resolution BSDF methods */
-SDFunc			SDhandleMtx = {
+const SDFunc		SDhandleMtx = {
 				&SDgetMtxBSDF,
 				&SDqueryMtxProjSA,
 				&SDgetMtxCDist,