--- ray/src/cv/bsdfrbf.c	2013/10/21 21:48:42	2.13
+++ ray/src/cv/bsdfrbf.c	2019/04/09 22:39:33	2.32
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfrbf.c,v 2.13 2013/10/21 21:48:42 greg Exp $";
+static const char RCSid[] = "$Id: bsdfrbf.c,v 2.32 2019/04/09 22:39:33 greg Exp $";
 #endif
 /*
  * Radial basis function representation for BSDF data.
@@ -29,21 +29,68 @@ static const char RCSid[] = "$Id: bsdfrbf.c,v 2.13 201
 #include "bsdfrep.h"
 
 #ifndef RSCA
-#define RSCA		2.2		/* radius scaling factor (empirical) */
+#define RSCA		2.0		/* radius scaling factor (empirical) */
 #endif
+#ifndef MAXSLOPE
+#define MAXSLOPE	200.0		/* maximum slope for smooth region */
+#endif
 #ifndef SMOOTH_MSE
 #define SMOOTH_MSE	5e-5		/* acceptable mean squared error */
 #endif
 #ifndef SMOOTH_MSER
-#define SMOOTH_MSER	0.07		/* acceptable relative MSE */
+#define SMOOTH_MSER	0.0025		/* acceptable relative MSE */
 #endif
 #define MAX_RAD		(GRIDRES/8)	/* maximum lobe radius */
 
 #define	RBFALLOCB	10		/* RBF allocation block size */
 
-				/* our loaded grid for this incident angle */
+					/* loaded grid or comparison DSFs */
 GRIDVAL			dsf_grid[GRIDRES][GRIDRES];
+					/* allocated chrominance sums if any */
+float			(*spec_grid)[GRIDRES][GRIDRES];
+int			nspec_grid = 0;
 
+/* Set up visible spectrum sampling */
+void
+set_spectral_samples(int nspec)
+{
+	if (rbf_colorimetry == RBCunknown) {
+		if (nspec_grid > 0) {
+			free(spec_grid);
+			spec_grid = NULL;
+			nspec_grid = 0;
+		}
+		if (nspec == 1) {
+			rbf_colorimetry = RBCphotopic;
+			return;
+		}
+		if (nspec == 3) {
+			rbf_colorimetry = RBCtristimulus;
+			spec_grid = (float (*)[GRIDRES][GRIDRES])calloc(
+					2*GRIDRES*GRIDRES, sizeof(float) );
+			if (spec_grid == NULL)
+				goto mem_error;
+			nspec_grid = 2;
+			return;
+		}
+		fprintf(stderr,
+			"%s: only 1 or 3 spectral samples currently supported\n",
+				progname);
+		exit(1);
+	}
+	if (nspec != nspec_grid+1) {
+		fprintf(stderr,
+			"%s: number of spectral samples cannot be changed\n",
+				progname);
+		exit(1);
+	}
+	return;
+mem_error:
+	fprintf(stderr, "%s: out of memory in set_spectral_samples()\n",
+			progname);
+	exit(1);
+}
+
 /* Start new DSF input grid */
 void
 new_bsdf_data(double new_theta, double new_phi)
@@ -51,48 +98,59 @@ new_bsdf_data(double new_theta, double new_phi)
 	if (!new_input_direction(new_theta, new_phi))
 		exit(1);
 	memset(dsf_grid, 0, sizeof(dsf_grid));
+	if (nspec_grid > 0)
+		memset(spec_grid, 0, sizeof(float)*GRIDRES*GRIDRES*nspec_grid);
 }
 
 /* Add BSDF data point */
 void
-add_bsdf_data(double theta_out, double phi_out, double val, int isDSF)
+add_bsdf_data(double theta_out, double phi_out, const double val[], int isDSF)
 {
 	FVECT	ovec;
+	double	cfact, Yval;
 	int	pos[2];
 
+	if (nspec_grid > 2) {
+		fprintf(stderr, "%s: unsupported color space\n", progname);
+		exit(1);
+	}
 	if (!output_orient)		/* check output orientation */
 		output_orient = 1 - 2*(theta_out > 90.);
 	else if (output_orient > 0 ^ theta_out < 90.) {
-		fputs("Cannot handle output angles on both sides of surface\n",
-				stderr);
+		fprintf(stderr,
+		"%s: cannot handle output angles on both sides of surface\n",
+				progname);
 		exit(1);
 	}
 	ovec[2] = sin((M_PI/180.)*theta_out);
 	ovec[0] = cos((M_PI/180.)*phi_out) * ovec[2];
 	ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
 	ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
+					/* BSDF to DSF correction */
+	cfact = isDSF ? 1. : ovec[2];
 
-	if (val <= 0)			/* truncate to zero */
-		val = 0;
-	else if (!isDSF)
-		val *= ovec[2];		/* convert from BSDF to DSF */
-
+	Yval = cfact * val[rbf_colorimetry==RBCtristimulus];
 					/* update BSDF histogram */
-	if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
-		++bsdf_hist[histndx(val/ovec[2])];
+	if (BSDF2SML*ovec[2] < Yval && Yval < BSDF2BIG*ovec[2])
+		++bsdf_hist[histndx(Yval/ovec[2])];
 
 	pos_from_vec(pos, ovec);
 
-	dsf_grid[pos[0]][pos[1]].vsum += val;
-	dsf_grid[pos[0]][pos[1]].nval++;
+	dsf_grid[pos[0]][pos[1]].sum.v += Yval;
+	dsf_grid[pos[0]][pos[1]].sum.n++;
+					/* add in X and Z values */
+	if (rbf_colorimetry == RBCtristimulus) {
+		spec_grid[0][pos[0]][pos[1]] += cfact * val[0];
+		spec_grid[1][pos[0]][pos[1]] += cfact * val[2];
+	}
 }
 
 /* Compute minimum BSDF from histogram (does not clear) */
 static void
 comp_bsdf_min()
 {
-	int	cnt;
-	int	i, target;
+	unsigned long	cnt, target;
+	int		i;
 
 	cnt = 0;
 	for (i = HISTLEN; i--; )
@@ -116,7 +174,7 @@ empty_region(int x0, int x1, int y0, int y1)
 
 	for (x = x0; x < x1; x++)
 	    for (y = y0; y < y1; y++)
-		if (dsf_grid[x][y].nval)
+		if (dsf_grid[x][y].sum.n)
 			return(0);
 	return(1);
 }
@@ -134,8 +192,8 @@ smooth_region(int x0, int x1, int y0, int y1)
 	memset(xvec, 0, sizeof(xvec));
 	for (x = x0; x < x1; x++)
 	    for (y = y0; y < y1; y++)
-		if ((n = dsf_grid[x][y].nval) > 0) {
-			double	z = dsf_grid[x][y].vsum;
+		if ((n = dsf_grid[x][y].sum.n) > 0) {
+			double	z = dsf_grid[x][y].sum.v;
 			rMtx[0][0] += x*x*(double)n;
 			rMtx[0][1] += x*y*(double)n;
 			rMtx[0][2] += x*(double)n;
@@ -147,18 +205,21 @@ smooth_region(int x0, int x1, int y0, int y1)
 			xvec[2] += z;
 		}
 	rMtx[1][0] = rMtx[0][1];
+	rMtx[2][0] = rMtx[0][2];
 	rMtx[2][1] = rMtx[1][2];
 	nvs = rMtx[2][2];
 	if (SDinvXform(rMtx, rMtx) != SDEnone)
-		return(0);
+		return(1);		/* colinear values */
 	A = DOT(rMtx[0], xvec);
 	B = DOT(rMtx[1], xvec);
+	if (A*A + B*B > MAXSLOPE*MAXSLOPE)	/* too steep? */
+		return(0);
 	C = DOT(rMtx[2], xvec);
 	sqerr = 0.0;			/* compute mean squared error */
 	for (x = x0; x < x1; x++)
 	    for (y = y0; y < y1; y++)
-		if ((n = dsf_grid[x][y].nval) > 0) {
-			double	d = A*x + B*y + C - dsf_grid[x][y].vsum/n;
+		if ((n = dsf_grid[x][y].sum.n) > 0) {
+			double	d = A*x + B*y + C - dsf_grid[x][y].sum.v/n;
 			sqerr += n*d*d;
 		}
 	if (sqerr <= nvs*SMOOTH_MSE)	/* below absolute MSE threshold? */
@@ -168,31 +229,58 @@ smooth_region(int x0, int x1, int y0, int y1)
 }
 
 /* Create new lobe based on integrated samples in region */
-static void
+static int
 create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
 {
 	double	vtot = 0.0;
+	double	CIEXtot = 0.0, CIEZtot = 0.0;
 	int	nv = 0;
+	double	wxsum = 0.0, wysum = 0.0, wtsum = 0.0;
 	double	rad;
 	int	x, y;
 					/* compute average for region */
 	for (x = x0; x < x1; x++)
-	    for (y = y0; y < y1; y++) {
-		vtot += dsf_grid[x][y].vsum;
-		nv += dsf_grid[x][y].nval;
-	    }
+	    for (y = y0; y < y1; y++)
+		if (dsf_grid[x][y].sum.n) {
+		    const double	v = dsf_grid[x][y].sum.v;
+		    const int		n = dsf_grid[x][y].sum.n;
+
+		    if (v > 0) {
+			const double	wt = v / (double)n;
+			wxsum += wt * x;
+			wysum += wt * y;
+			wtsum += wt;
+		    }
+		    vtot += v;
+		    nv += n;
+		    if (rbf_colorimetry == RBCtristimulus) {
+			CIEXtot += spec_grid[0][x][y];
+			CIEZtot += spec_grid[1][x][y];
+		    }
+		}
 	if (!nv) {
 		fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
 				progname);
 		exit(1);
 	}
+	if (vtot <= 0)			/* only create positive lobes */
+		return(0);
+					/* assign color */
+	if (rbf_colorimetry == RBCtristimulus) {
+		const double	df = 1.0 / (CIEXtot + vtot + CIEZtot);
+		C_COLOR		cclr;
+		c_cset(&cclr, CIEXtot*df, vtot*df);
+		rvp->chroma = c_encodeChroma(&cclr);
+	} else
+		rvp->chroma = c_dfchroma;
 					/* peak value based on integral */
 	vtot *= (x1-x0)*(y1-y0)*(2.*M_PI/GRIDRES/GRIDRES)/(double)nv;
 	rad = (RSCA/(double)GRIDRES)*(x1-x0);
 	rvp->peak =  vtot / ((2.*M_PI) * rad*rad);
-	rvp->crad = ANG2R(rad);
-	rvp->gx = (x0+x1)>>1;
-	rvp->gy = (y0+y1)>>1;
+	rvp->crad = ANG2R(rad);		/* put peak at centroid */
+	rvp->gx = (int)(wxsum/wtsum + .5);
+	rvp->gy = (int)(wysum/wtsum + .5);
+	return(1);
 }
 
 /* Recursive function to build radial basis function representation */
@@ -223,7 +311,8 @@ build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, in
 	if (!nleaves)			/* nothing but branches? */
 		return(nadded);
 					/* combine 4 leaves into 1? */
-	if (nleaves == 4 && x1-x0 < MAX_RAD && smooth_region(x0, x1, y0, y1))
+	if ((nleaves == 4) & (x1-x0 <= MAX_RAD) &&
+			smooth_region(x0, x1, y0, y1))
 		return(0);
 					/* need more array space? */
 	if ((*np+nleaves-1)>>RBFALLOCB != (*np-1)>>RBFALLOCB) {
@@ -233,15 +322,18 @@ build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, in
 			return(-1);
 	}
 					/* create lobes for leaves */
-	if (!branched[0])
-		create_lobe(*arp+(*np)++, x0, xmid, y0, ymid);
-	if (!branched[1])
-		create_lobe(*arp+(*np)++, xmid, x1, y0, ymid);
-	if (!branched[2])
-		create_lobe(*arp+(*np)++, x0, xmid, ymid, y1);
-	if (!branched[3])
-		create_lobe(*arp+(*np)++, xmid, x1, ymid, y1);
-	nadded += nleaves;
+	if (!branched[0] && create_lobe(*arp+*np, x0, xmid, y0, ymid)) {
+		++(*np); ++nadded;
+	}
+	if (!branched[1] && create_lobe(*arp+*np, xmid, x1, y0, ymid)) {
+		++(*np); ++nadded;
+	}
+	if (!branched[2] && create_lobe(*arp+*np, x0, xmid, ymid, y1)) {
+		++(*np); ++nadded;
+	}
+	if (!branched[3] && create_lobe(*arp+*np, xmid, x1, ymid, y1)) {
+		++(*np); ++nadded;
+	}
 	return(nadded);
 }
 
@@ -256,8 +348,14 @@ make_rbfrep()
 	comp_bsdf_min();
 				/* create RBF node list */
 	rbfarr = NULL; nn = 0;
-	if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0)
-		goto memerr;
+	if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0) {
+		if (nn)
+			goto memerr;
+		fprintf(stderr,
+			"%s: warning - skipping bad incidence (%.1f,%.1f)\n",
+				progname, theta_in_deg, phi_in_deg);
+		return(NULL);
+	}
 				/* (re)allocate RBF array */
 	newnode = (RBFNODE *)realloc(rbfarr,
 			sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
@@ -284,7 +382,6 @@ make_rbfrep()
 			newnode->vtotal);
 #endif
 	insert_dsf(newnode);
-
 	return(newnode);
 memerr:
 	fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);