--- ray/src/common/interp2d.c	2013/02/09 17:39:21	2.2
+++ ray/src/common/interp2d.c	2013/02/12 02:56:05	2.6
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: interp2d.c,v 2.2 2013/02/09 17:39:21 greg Exp $";
+static const char RCSid[] = "$Id: interp2d.c,v 2.6 2013/02/12 02:56:05 greg Exp $";
 #endif
 /*
  * General interpolation method for unstructured values on 2-D plane.
@@ -9,29 +9,30 @@ static const char RCSid[] = "$Id: interp2d.c,v 2.2 201
 
 #include "copyright.h"
 
-/*************************************************************
+/***************************************************************
  * This is a general method for 2-D interpolation similar to
  * radial basis functions but allowing for a good deal of local
  * anisotropy in the point distribution.  Each sample point
  * is examined to determine the closest neighboring samples in
  * each of NI2DIR surrounding directions.  To speed this
- * calculation, we sort the data into 3 half-planes and
- * perform simple tests to see which neighbor is closest in
- * a each direction.  Once we have our approximate neighborhood
- * for a sample, we can use it in a Gaussian weighting scheme
- * with anisotropic surround.  This gives us a fairly smooth
- * interpolation however the sample points may be initially
- * distributed.  Evaluation is accelerated by use of a fast
- * approximation to the atan2(y,x) function.
- **************************************************************/
+ * calculation, we sort the data into half-planes and apply
+ * simple tests to see which neighbor is closest in each
+ * direction.  Once we have our approximate neighborhood
+ * for a sample, we can use it in a modified Gaussian weighting
+ * with allowing local anisotropy.  Harmonic weighting is added
+ * 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.
+ ****************************************************************/
 
 #include <stdio.h>
 #include <stdlib.h>
 #include "rtmath.h"
 #include "interp2d.h"
 
-#define DECODE_RAD(ip,er)	((ip)->rmin*(1. + .5*(er)))
-#define ENCODE_RAD(ip,r)	((int)(2.*(r)/(ip)->rmin) - 2)
+#define DECODE_DIA(ip,ed)	((ip)->dmin*(1. + .5*(ed)))
+#define ENCODE_DIA(ip,d)	((int)(2.*(d)/(ip)->dmin) - 2)
 
 /* Sample order (private) */
 typedef struct {
@@ -53,9 +54,9 @@ interp2_alloc(int nsamps)
 		return(NULL);
 
 	nip->ns = nsamps;
-	nip->rmin = .5;		/* default radius minimum */
+	nip->dmin = 1;		/* default minimum diameter */
 	nip->smf = NI2DSMF;	/* default smoothing factor */
-	nip->ra = NULL;
+	nip->da = NULL;
 				/* caller must assign spt[] array */
 	return(nip);
 }
@@ -72,9 +73,9 @@ interp2_realloc(INTERP2 *ip, int nsamps)
 	}
 	if (nsamps == ip->ns);
 		return(ip);
-	if (ip->ra != NULL) {	/* will need to recompute distribution */
-		free(ip->ra);
-		ip->ra = NULL;
+	if (ip->da != NULL) {	/* will need to recompute distribution */
+		free(ip->da);
+		ip->da = NULL;
 	}
 	ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
 	if (ip == NULL)
@@ -83,6 +84,29 @@ interp2_realloc(INTERP2 *ip, int nsamps)
 	return(ip);
 }
 
+/* Set minimum distance under which samples will start to merge */
+void
+interp2_spacing(INTERP2 *ip, double mind)
+{
+	if (mind <= 0)
+		return;
+	if ((.998*ip->dmin <= mind) && (mind <= 1.002*ip->dmin))
+		return;
+	if (ip->da != NULL) {	/* will need to recompute distribution */
+		free(ip->da);
+		ip->da = NULL;
+	}
+	ip->dmin = mind;
+}
+
+/* Modify smoothing parameter by the given factor */
+void
+interp2_smooth(INTERP2 *ip, double sf)
+{
+	if ((ip->smf *= sf) < NI2DSMF)
+		ip->smf = NI2DSMF;
+}
+
 /* private call-back to sort position index */
 static int
 cmp_spos(const void *p1, const void *p2)
@@ -112,17 +136,17 @@ sort_samples(SAMPORD *sord, const INTERP2 *ip, double 
 	qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
 }
 
-/* private routine to encode radius with range checks */
+/* private routine to encode sample diameter with range checks */
 static int
-encode_radius(const INTERP2 *ip, double r)
+encode_diameter(const INTERP2 *ip, double d)
 {
-	const int	er = ENCODE_RAD(ip, r);
+	const int	ed = ENCODE_DIA(ip, d);
 
-	if (er <= 0)
+	if (ed <= 0)
 		return(0);
-	if (er >= 0xffff)
+	if (ed >= 0xffff)
 		return(0xffff);
-	return(er);
+	return(ed);
 }
 
 /* (Re)compute anisotropic basis function interpolant (normally automatic) */
@@ -133,13 +157,13 @@ interp2_analyze(INTERP2 *ip)
 	int	*rightrndx, *leftrndx, *endrndx;
 	int	bd;
 					/* sanity checks */
-	if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
+	if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
 		return(0);
 					/* need to allocate? */
-	if (ip->ra == NULL) {
-		ip->ra = (unsigned short (*)[NI2DIR])malloc(
+	if (ip->da == NULL) {
+		ip->da = (unsigned short (*)[NI2DIR])malloc(
 				sizeof(unsigned short)*NI2DIR*ip->ns);
-		if (ip->ra == NULL)
+		if (ip->da == NULL)
 			return(0);
 	}
 					/* get temporary arrays */
@@ -182,24 +206,23 @@ interp2_analyze(INTERP2 *ip)
 	    sort_samples(sortord, ip, ang);
 					/* find nearest neighbors each side */
 	    for (i = ip->ns; i--; ) {
-		const int	rpos = rightrndx[sortord[i].si];
-		const int	lpos = leftrndx[sortord[i].si];
+		const int	ii = sortord[i].si;
 		int		j;
-					/* preload with large radius */
-		ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
-			    .25*(sortord[ip->ns-1].dm - sortord[0].dm));
+					/* 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));
 		for (j = i; ++j < ip->ns; )	/* nearest above */
-		    if (rightrndx[sortord[j].si] > rpos &&
-				    leftrndx[sortord[j].si] < lpos) {
-			ip->ra[i][bd] = encode_radius(ip,
-					.5*(sortord[j].dm - sortord[i].dm));
+		    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
+				    leftrndx[sortord[j].si] < leftrndx[ii]) {
+			ip->da[ii][bd] = encode_diameter(ip,
+						sortord[j].dm - sortord[i].dm);
 			break;
 		    }
 		for (j = i; j-- > 0; )		/* nearest below */
-		    if (rightrndx[sortord[j].si] < rpos &&
-				    leftrndx[sortord[j].si] > lpos) {
-			ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
-					.5*(sortord[i].dm - sortord[j].dm));
+		    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
+				    leftrndx[sortord[j].si] > leftrndx[ii]) {
+			ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
+						sortord[i].dm - sortord[j].dm);
 			break;
 		    }
 	    }
@@ -211,11 +234,11 @@ interp2_analyze(INTERP2 *ip)
 	return(1);
 }
 
-/* private call returns log of raw weight for a particular sample */
+/* private call returns raw weight for a particular sample */
 static double
-get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
+get_wt(const INTERP2 *ip, const int i, double x, double y)
 {
-	double	dir, rd;
+	double	dir, rd, r2, d2;
 	int	ri;
 				/* get relative direction */
 	x -= ip->spt[i][0];
@@ -226,10 +249,14 @@ get_ln_wt(const INTERP2 *ip, const int i, double x, do
 	rd = dir * (NI2DIR/2./PI);
 	ri = (int)rd;
 	rd -= (double)ri;
-	rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
-	rd = ip->smf * DECODE_RAD(ip, rd);
-				/* return log of Gaussian weight */
-	return( (x*x + y*y) / (-2.*rd*rd) );
+	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 */
@@ -242,17 +269,12 @@ interp2_weights(float wtv[], INTERP2 *ip, double x, do
 	if ((wtv == NULL) | (ip == NULL))
 		return(0);
 					/* need to compute interpolant? */
-	if (ip->ra == NULL && !interp2_analyze(ip))
+	if (ip->da == NULL && !interp2_analyze(ip))
 		return(0);
 
 	wnorm = 0;			/* compute raw weights */
 	for (i = ip->ns; i--; ) {
-		double	wt = get_ln_wt(ip, i, x, y);
-		if (wt < -21.) {
-			wtv[i] = 0;	/* ignore weights < 1e-9 */
-			continue;
-		}
-		wt = exp(wt);		/* Gaussian weight */
+		double	wt = get_wt(ip, i, x, y);
 		wtv[i] = wt;
 		wnorm += wt;
 	}
@@ -276,13 +298,13 @@ interp2_topsamp(float wt[], int si[], const int n, INT
 	if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
 		return(0);
 					/* need to compute interpolant? */
-	if (ip->ra == NULL && !interp2_analyze(ip))
+	if (ip->da == NULL && !interp2_analyze(ip))
 		return(0);
 					/* identify top n weights */
 	for (i = ip->ns; i--; ) {
-		const double	lnwt = get_ln_wt(ip, i, x, y);
+		const double	wti = get_wt(ip, i, x, y);
 		for (j = nn; j > 0; j--) {
-			if (wt[j-1] >= lnwt)
+			if (wt[j-1] >= wti)
 				break;
 			if (j < n) {
 				wt[j] = wt[j-1];
@@ -290,17 +312,14 @@ interp2_topsamp(float wt[], int si[], const int n, INT
 			}
 		}
 		if (j < n) {		/* add/insert sample */
-			wt[j] = lnwt;
+			wt[j] = wti;
 			si[j] = i;
 			nn += (nn < n);
 		}
 	}
-	wnorm = 0;			/* exponentiate and normalize */
-	for (j = nn; j--; ) {
-		double	dwt = exp(wt[j]);
-		wt[j] = dwt;
-		wnorm += dwt;
-	}
+	wnorm = 0;			/* normalize sample weights */
+	for (j = nn; j--; )
+		wnorm += wt[j];
 	if (wnorm <= 0)
 		return(0);
 	wnorm = 1./wnorm;
@@ -315,7 +334,7 @@ interp2_free(INTERP2 *ip)
 {
 	if (ip == NULL)
 		return;
-	if (ip->ra != NULL)
-		free(ip->ra);
+	if (ip->da != NULL)
+		free(ip->da);
 	free(ip);
 }