--- ray/src/common/interp2d.c	2013/02/12 18:13:28	2.8
+++ ray/src/common/interp2d.c	2013/02/15 01:26:47	2.11
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: interp2d.c,v 2.8 2013/02/12 18:13:28 greg Exp $";
+static const char RCSid[] = "$Id: interp2d.c,v 2.11 2013/02/15 01:26:47 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.8 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 an array
+ * of flags indicating where weights are (nearly) zero.
  ****************************************************************/
 
 #include <stdio.h>
@@ -155,17 +156,39 @@ 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;
 	}
+	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 *)calloc( ip->ns,
+					sizeof(struct interp2_samp) );
+	if (ip->da == NULL)
+		return(0);
 					/* get temporary arrays */
 	sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
 	rightrndx = (int *)malloc(sizeof(int)*ip->ns);
@@ -178,7 +201,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,19 +231,20 @@ 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;
 		    }
@@ -234,12 +257,12 @@ interp2_analyze(INTERP2 *ip)
 	return(1);
 }
 
-/* private call returns raw weight for a particular sample */
-static double
-get_wt(const INTERP2 *ip, const int i, double x, double y)
+/* 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;
+	double		dir, rd, r2, d2;
+	int		ri;
 				/* get relative direction */
 	x -= ip->spt[i][0];
 	y -= ip->spt[i][1];
@@ -249,7 +272,7 @@ get_wt(const INTERP2 *ip, const int i, double x, doubl
 	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 = (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;
@@ -259,25 +282,78 @@ get_wt(const INTERP2 *ip, const int i, double x, doubl
 	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	ingrid = 1;
+
+	if (ip == NULL)		/* paranoia */
+		return(-1);
+				/* need to compute interpolant? */
+	if (ip->da == NULL && !interp2_analyze(ip))
+		return(-1);
+				/* get grid position */
+	fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
+	if (fgi[0] >= NI2DIM) {
+		fgi[0] = NI2DIM-1;
+		ingrid = 0;
+	} else if (fgi[0] < 0) {
+		fgi[0] = 0;
+		ingrid = 0;
+	}
+	fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
+	if (fgi[1] >= NI2DIM) {
+		fgi[1] = NI2DIM-1;
+		ingrid = 0;
+	} else if (fgi[1] < 0) {
+		fgi[1] = 0;
+		ingrid = 0;
+	}
+	return(ingrid);
+}
+
+/* private call to set black flag if not too close to the given sample */
+static void
+setblk(INTERP2 *ip, const int i, const int gi[2])
+{
+	double	dx = (gi[0]+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
+					ip->smin[0] - ip->spt[i][0];
+	double	dy = (gi[1]+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
+					ip->smin[1] - ip->spt[i][1];
+
+	if (dx*dx + dy*dy > 2.*ip->grid2)
+		ip->da[i].blkflg[gi[1]] |= 1<<gi[0];
+}
+
+#define chkblk(ip,i,gi)	((ip)->da[i].blkflg[(gi)[1]]>>(gi)[0] & 1)
+
 /* 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	ingrid;
 	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 ((ingrid = 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 (chkblk(ip, i, fgi)) {
+		wtv[i] = 0;
+	    } else {
+		double	wt = interp2_wti(ip, i, x, y);
 		wtv[i] = wt;
 		wnorm += wt;
-	}
+		if (wt <= 1e-9 && ingrid)
+			setblk(ip, i, fgi);
+	    }
 	if (wnorm <= 0)			/* too far from all our samples! */
 		return(0);
 	wnorm = 1./wnorm;		/* normalize weights */
@@ -292,17 +368,25 @@ int
 interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
 {
 	int	nn = 0;
+	int	fgi[2];
+	int	ingrid;
 	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 ((ingrid = 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);
+	for (i = ip->ns; i--; )
+	    if (!chkblk(ip, i, fgi)) {
+		const double	wti = interp2_wti(ip, i, x, y);
+		if (wti <= 1e-9) {
+			if (ingrid)
+				setblk(ip, i, fgi);
+			continue;
+		}
 		for (j = nn; j > 0; j--) {
 			if (wt[j-1] >= wti)
 				break;
@@ -316,7 +400,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];