--- ray/src/common/interp2d.c	2013/02/12 18:41:39	2.9
+++ ray/src/common/interp2d.c	2014/06/06 00:56:42	2.14
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: interp2d.c,v 2.9 2013/02/12 18:41:39 greg Exp $";
+static const char RCSid[] = "$Id: interp2d.c,v 2.14 2014/06/06 00:56:42 greg Exp $";
 #endif
 /*
  * General interpolation method for unstructured values on 2-D plane.
@@ -23,7 +23,8 @@ static const char RCSid[] = "$Id: interp2d.c,v 2.9 201
  * to reduce the influence of distant neighbors.  This yields a
  * smooth interpolation regardless of how the sample points are
  * initially distributed.  Evaluation is accelerated by use of
- * a fast approximation to the atan2(y,x) function.
+ * a fast approximation to the atan2(y,x) function and a low-res
+ * map indicating where sample weights are significant.
  ****************************************************************/
 
 #include <stdio.h>
@@ -40,6 +41,19 @@ typedef struct {
 	float	dm;		/* distance measure in this direction */
 } SAMPORD;
 
+/* private routine to encode sample diameter with range checks */
+static int
+encode_diameter(const INTERP2 *ip, double d)
+{
+	const int	ed = ENCODE_DIA(ip, d);
+
+	if (ed <= 0)
+		return(0);
+	if (ed >= 0xffff)
+		return(0xffff);
+	return(ed);
+}
+
 /* Allocate a new set of interpolation samples (caller assigns spt[] array) */
 INTERP2	*
 interp2_alloc(int nsamps)
@@ -99,12 +113,131 @@ interp2_spacing(INTERP2 *ip, double mind)
 	ip->dmin = mind;
 }
 
+/* Compute unnormalized weight for a position relative to a sample */
+double
+interp2_wti(INTERP2 *ip, const int i, double x, double y)
+{
+	double		dir, rd, r2, d2;
+	int		ri;
+				/* get relative direction */
+	x -= ip->spt[i][0];
+	y -= ip->spt[i][1];
+	dir = atan2a(y, x);
+	dir += 2.*PI*(dir < 0);
+				/* linear radius interpolation */
+	rd = dir * (NI2DIR/2./PI);
+	ri = (int)rd;
+	rd -= (double)ri;
+	rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
+	rd = ip->smf * DECODE_DIA(ip, rd);
+	r2 = 2.*rd*rd;
+	d2 = x*x + y*y;
+	if (d2 > 21.*r2)	/* result would be < 1e-9 */
+		return(.0);
+				/* Gaussian times harmonic weighting */
+	return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
+}
+
+/* private call to get grid flag index */
+static int
+interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
+{
+	int	inbounds = 0;
+
+	if (ip == NULL)		/* paranoia */
+		return(-1);
+				/* need to compute interpolant? */
+	if (ip->da == NULL && !interp2_analyze(ip))
+		return(-1);
+				/* get x & y grid positions */
+	fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
+
+	if (fgi[0] >= NI2DIM)
+		fgi[0] = NI2DIM-1;
+	else if (fgi[0] < 0)
+		fgi[0] = 0;
+	else
+		++inbounds;
+
+	fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
+
+	if (fgi[1] >= NI2DIM)
+		fgi[1] = NI2DIM-1;
+	else if (fgi[1] < 0)
+		fgi[1] = 0;
+	else
+		++inbounds;
+
+	return(inbounds == 2);
+}
+
+#define setflg(fl,xi,yi)	((fl)[yi] |= 1<<(xi))
+
+#define chkflg(fl,xi,yi)	((fl)[yi]>>(xi) & 1)
+
+/* private flood function to determine sample influence */
+static void
+influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
+			int xfi, int yfi)
+{
+	double	gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
+					ip->smin[0];
+	double	gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
+					ip->smin[1];
+	double	dx = gx - ip->spt[i][0];
+	double	dy = gy - ip->spt[i][1];
+
+	setflg(visited, xfi, yfi);
+
+	if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
+		return;
+
+	setflg(ip->da[i].infl, xfi, yfi);
+
+	if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
+		influence_flood(ip, i, visited, xfi-1, yfi);
+
+	if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
+		influence_flood(ip, i, visited, xfi, yfi-1);
+
+	if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
+		influence_flood(ip, i, visited, xfi+1, yfi);
+
+	if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
+		influence_flood(ip, i, visited, xfi, yfi+1);
+}
+
+/* private call to compute sample influence maps */
+static void
+map_influence(INTERP2 *ip)
+{
+	unsigned short	visited[NI2DIM];
+	int		fgi[2];
+	int		i, j;
+
+	for (i = ip->ns; i--; ) {
+		for (j = NI2DIM; j--; ) {
+			ip->da[i].infl[j] = 0;
+			visited[j] = 0;
+		}
+		interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
+
+		influence_flood(ip, i, visited, fgi[0], fgi[1]);
+	}
+}
+
 /* Modify smoothing parameter by the given factor */
 void
 interp2_smooth(INTERP2 *ip, double sf)
 {
+	float	old_smf = ip->smf;
+
 	if ((ip->smf *= sf) < NI2DSMF)
 		ip->smf = NI2DSMF;
+					/* need to recompute influence maps? */
+	if (ip->da != NULL && (old_smf*.85 > ip->smf) |
+				(ip->smf > old_smf*1.15))
+		map_influence(ip);
 }
 
 /* private call-back to sort position index */
@@ -133,40 +266,45 @@ sort_samples(SAMPORD *sord, const INTERP2 *ip, double 
 		sord[i].si = i;
 		sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
 	}
-	qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
+	qsort(sord, ip->ns, sizeof(SAMPORD), cmp_spos);
 }
 
-/* private routine to encode sample diameter with range checks */
-static int
-encode_diameter(const INTERP2 *ip, double d)
-{
-	const int	ed = ENCODE_DIA(ip, d);
-
-	if (ed <= 0)
-		return(0);
-	if (ed >= 0xffff)
-		return(0xffff);
-	return(ed);
-}
-
 /* (Re)compute anisotropic basis function interpolant (normally automatic) */
 int
 interp2_analyze(INTERP2 *ip)
 {
 	SAMPORD	*sortord;
 	int	*rightrndx, *leftrndx, *endrndx;
-	int	bd;
+	int	i, bd;
 					/* sanity checks */
-	if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
+	if (ip == NULL)
 		return(0);
-					/* need to allocate? */
-	if (ip->da == NULL) {
-		ip->da = (unsigned short (*)[NI2DIR])malloc(
-				sizeof(unsigned short)*NI2DIR*ip->ns);
-		if (ip->da == NULL)
-			return(0);
+	if (ip->da != NULL) {		/* free previous data if any */
+		free(ip->da);
+		ip->da = NULL;
 	}
-					/* get temporary arrays */
+	if ((ip->ns <= 1) | (ip->dmin <= 0))
+		return(0);
+					/* compute sample domain */
+	ip->smin[0] = ip->smin[1] = FHUGE;
+	ip->smax[0] = ip->smax[1] = -FHUGE;
+	for (i = ip->ns; i--; ) {
+		if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
+		if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
+		if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
+		if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
+	}
+	ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
+			(ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
+				(1./NI2DIM/NI2DIM);
+	if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
+		return(0);
+					/* allocate analysis data */
+	ip->da = (struct interp2_samp *)malloc(
+				sizeof(struct interp2_samp)*ip->ns );
+	if (ip->da == NULL)
+		return(0);
+					/* allocate temporary arrays */
 	sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
 	rightrndx = (int *)malloc(sizeof(int)*ip->ns);
 	leftrndx = (int *)malloc(sizeof(int)*ip->ns);
@@ -178,7 +316,6 @@ interp2_analyze(INTERP2 *ip)
 	for (bd = 0; bd < NI2DIR/2; bd++) {
 	    const double	ang = 2.*PI/NI2DIR*bd;
 	    int			*sptr;
-	    int			i;
 					/* create right reverse index */
 	    if (bd) {			/* re-use from previous iteration? */
 		sptr = rightrndx;
@@ -209,75 +346,56 @@ interp2_analyze(INTERP2 *ip)
 		const int	ii = sortord[i].si;
 		int		j;
 					/* preload with large radii */
-		ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
-			    .5*(sortord[ip->ns-1].dm - sortord[0].dm));
+		ip->da[ii].dia[bd] =
+		ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
+				.5*(sortord[ip->ns-1].dm - sortord[0].dm));
 		for (j = i; ++j < ip->ns; )	/* nearest above */
 		    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
 				    leftrndx[sortord[j].si] < leftrndx[ii]) {
-			ip->da[ii][bd] = encode_diameter(ip,
+			ip->da[ii].dia[bd] = encode_diameter(ip,
 						sortord[j].dm - sortord[i].dm);
 			break;
 		    }
 		for (j = i; j-- > 0; )		/* nearest below */
 		    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
 				    leftrndx[sortord[j].si] > leftrndx[ii]) {
-			ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
+			ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
 						sortord[i].dm - sortord[j].dm);
 			break;
 		    }
 	    }
 	}
-	free(sortord);			/* clean up */
+	free(sortord);			/* release temp arrays */
 	free(rightrndx);
 	free(leftrndx);
 	free(endrndx);
-	return(1);
+					/* map sample influence areas */
+	map_influence(ip);
+	return(1);			/* all done */
 }
 
-/* private call returns raw weight for a particular sample */
-static double
-get_wt(const INTERP2 *ip, const int i, double x, double y)
-{
-	double	dir, rd, r2, d2;
-	int	ri;
-				/* get relative direction */
-	x -= ip->spt[i][0];
-	y -= ip->spt[i][1];
-	dir = atan2a(y, x);
-	dir += 2.*PI*(dir < 0);
-				/* linear radius interpolation */
-	rd = dir * (NI2DIR/2./PI);
-	ri = (int)rd;
-	rd -= (double)ri;
-	rd = (1.-rd)*ip->da[i][ri] + rd*ip->da[i][(ri+1)%NI2DIR];
-	rd = ip->smf * DECODE_DIA(ip, rd);
-	r2 = 2.*rd*rd;
-	d2 = x*x + y*y;
-	if (d2 > 21.*r2)	/* result would be < 1e-9 */
-		return(.0);
-				/* Gaussian times harmonic weighting */
-	return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
-}
-
 /* Assign full set of normalized weights to interpolate the given position */
 int
 interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
 {
 	double	wnorm;
+	int	fgi[2];
 	int	i;
 
-	if ((wtv == NULL) | (ip == NULL))
+	if (wtv == NULL)
 		return(0);
-					/* need to compute interpolant? */
-	if (ip->da == NULL && !interp2_analyze(ip))
+					/* get flag position */
+	if (interp2_flagpos(fgi, ip, x, y) < 0)
 		return(0);
 
 	wnorm = 0;			/* compute raw weights */
-	for (i = ip->ns; i--; ) {
-		double	wt = get_wt(ip, i, x, y);
+	for (i = ip->ns; i--; )
+	    if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
+		double	wt = interp2_wti(ip, i, x, y);
 		wtv[i] = wt;
 		wnorm += wt;
-	}
+	    } else
+		wtv[i] = 0;
 	if (wnorm <= 0)			/* too far from all our samples! */
 		return(0);
 	wnorm = 1./wnorm;		/* normalize weights */
@@ -292,19 +410,19 @@ int
 interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
 {
 	int	nn = 0;
+	int	fgi[2];
 	double	wnorm;
 	int	i, j;
 
-	if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
+	if ((n <= 0) | (wt == NULL) | (si == NULL))
 		return(0);
-					/* need to compute interpolant? */
-	if (ip->da == NULL && !interp2_analyze(ip))
+					/* get flag position */
+	if (interp2_flagpos(fgi, ip, x, y) < 0)
 		return(0);
 					/* identify top n weights */
-	for (i = ip->ns; i--; ) {
-		const double	wti = get_wt(ip, i, x, y);
-		if (wti <= 1e-9)
-			continue;
+	for (i = ip->ns; i--; )
+	    if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
+		const double	wti = interp2_wti(ip, i, x, y);
 		for (j = nn; j > 0; j--) {
 			if (wt[j-1] >= wti)
 				break;
@@ -318,7 +436,7 @@ interp2_topsamp(float wt[], int si[], const int n, INT
 			si[j] = i;
 			nn += (nn < n);
 		}
-	}
+	    }
 	wnorm = 0;			/* normalize sample weights */
 	for (j = nn; j--; )
 		wnorm += wt[j];