src/common/interp2d.c

#ifndef lint
static const char RCSid[] = "$Id: interp2d.c,v 2.4 2013/02/11 22:56:22 greg Exp $";
#endif
/*
 * General interpolation method for unstructured values on 2-D plane.
 *
 *      G.Ward Feb 2013
 */

#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 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_DIA(ip,ed)       ((ip)->dmin*(1. + .5*(ed)))
#define ENCODE_DIA(ip,d)        ((int)(2.*(d)/(ip)->dmin) - 2)

/* Sample order (private) */
typedef struct {
        int     si;             /* sample index */
        float   dm;             /* distance measure in this direction */
} SAMPORD;

/* Allocate a new set of interpolation samples (caller assigns spt[] array) */
INTERP2 *
interp2_alloc(int nsamps)
{
        INTERP2         *nip;

        if (nsamps <= 1)
                return(NULL);

        nip = (INTERP2 *)malloc(sizeof(INTERP2) + sizeof(float)*2*(nsamps-1));
        if (nip == NULL)
                return(NULL);

        nip->ns = nsamps;
        nip->dmin = 1;          /* default minimum diameter */
        nip->smf = NI2DSMF;     /* default smoothing factor */
        nip->da = NULL;
                                /* caller must assign spt[] array */
        return(nip);
}

/* Resize interpolation array (caller must assign any new values) */
INTERP2 *
interp2_realloc(INTERP2 *ip, int nsamps)
{
        if (ip == NULL)
                return(interp2_alloc(nsamps));
        if (nsamps <= 1) {
                interp2_free(ip);
                return(NULL);
        }
        if (nsamps == ip->ns);
                return(ip);
        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)
                return(NULL);
        ip->ns = 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)
{
        const SAMPORD   *so1 = (const SAMPORD *)p1;
        const SAMPORD   *so2 = (const SAMPORD *)p2;

        if (so1->dm > so2->dm)
                return 1;
        if (so1->dm < so2->dm)
                return -1;
        return 0;
}

/* private routine to order samples in a particular direction */
static void
sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang)
{
        const double    cosd = cos(ang);
        const double    sind = sin(ang);
        int             i;

        for (i = ip->ns; i--; ) {
                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);
}

/* 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;
                                        /* sanity checks */
        if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
                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);
        }
                                        /* get temporary arrays */
        sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
        rightrndx = (int *)malloc(sizeof(int)*ip->ns);
        leftrndx = (int *)malloc(sizeof(int)*ip->ns);
        endrndx = (int *)malloc(sizeof(int)*ip->ns);
        if ((sortord == NULL) | (rightrndx == NULL) |
                        (leftrndx == NULL) | (endrndx == NULL))
                return(0);
                                        /* run through bidirections */
        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;
                rightrndx = leftrndx;
                leftrndx = sptr;
            } else {                    /* else sort first half-plane */
                sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
                for (i = ip->ns; i--; )
                    rightrndx[sortord[i].si] = i;
                                        /* & store reverse order for later */
                for (i = ip->ns; i--; )
                    endrndx[sortord[i].si] = ip->ns-1 - i;
            }
                                        /* create new left reverse index */
            if (bd == NI2DIR/2 - 1) {   /* use order from first iteration? */
                sptr = leftrndx;
                leftrndx = endrndx;
                endrndx = sptr;
            } else {                    /* else compute new half-plane */
                sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
                for (i = ip->ns; i--; )
                    leftrndx[sortord[i].si] = i;
            }
                                        /* sort grid values in this direction */
            sort_samples(sortord, ip, ang);
                                        /* find nearest neighbors each side */
            for (i = ip->ns; i--; ) {
                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));
                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,
                                                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,
                                                sortord[i].dm - sortord[j].dm);
                        break;
                    }
            }
        }
        free(sortord);                  /* clean up */
        free(rightrndx);
        free(leftrndx);
        free(endrndx);
        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)
{
        double  dir, rd, 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);
        d2 = x*x + y*y;
                                /* Gaussian times harmonic weighting */
        return( exp(d2/(-2.*rd*rd)) * 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     i;

        if ((wtv == NULL) | (ip == NULL))
                return(0);
                                        /* need to compute interpolant? */
        if (ip->da == NULL && !interp2_analyze(ip))
                return(0);

        wnorm = 0;                      /* compute raw weights */
        for (i = ip->ns; i--; ) {
                double  wt = get_wt(ip, i, x, y);
                wtv[i] = wt;
                wnorm += wt;
        }
        if (wnorm <= 0)                 /* too far from all our samples! */
                return(0);
        wnorm = 1./wnorm;               /* normalize weights */
        for (i = ip->ns; i--; )
                wtv[i] *= wnorm;
        return(ip->ns);                 /* all done */
}


/* Get normalized weights and indexes for n best samples in descending order */
int
interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
{
        int     nn = 0;
        double  wnorm;
        int     i, j;

        if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
                return(0);
                                        /* need to compute interpolant? */
        if (ip->da == NULL && !interp2_analyze(ip))
                return(0);
                                        /* identify top n weights */
        for (i = ip->ns; i--; ) {
                const double    wti = get_wt(ip, i, x, y);
                for (j = nn; j > 0; j--) {
                        if (wt[j-1] >= wti)
                                break;
                        if (j < n) {
                                wt[j] = wt[j-1];
                                si[j] = si[j-1];
                        }
                }
                if (j < n) {            /* add/insert sample */
                        wt[j] = wti;
                        si[j] = i;
                        nn += (nn < n);
                }
        }
        wnorm = 0;                      /* normalize sample weights */
        for (j = nn; j--; )
                wnorm += wt[j];
        if (wnorm <= 0)
                return(0);
        wnorm = 1./wnorm;
        for (j = nn; j--; )
                wt[j] *= wnorm;
        return(nn);                     /* return actual # samples */
}

/* Free interpolant */
void
interp2_free(INTERP2 *ip)
{
        if (ip == NULL)
                return;
        if (ip->da != NULL)
                free(ip->da);
        free(ip);
}
Revision:	2.5
Committed:	Mon Feb 11 23:33:35 2013 UTC (11 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.4:	+24 -1 lines
Log Message:	Added accessor functions for minimum spacing and smoothing parameters
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.5	static const char RCSid[] = "$Id: interp2d.c,v 2.4 2013/02/11 22:56:22 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* General interpolation method for unstructured values on 2-D plane.
6			*
7			* G.Ward Feb 2013
8			*/
9
10			#include "copyright.h"
11
12	greg	2.4	/***************************************************************
13	greg	2.1	* This is a general method for 2-D interpolation similar to
14			* radial basis functions but allowing for a good deal of local
15			* anisotropy in the point distribution. Each sample point
16			* is examined to determine the closest neighboring samples in
17			* each of NI2DIR surrounding directions. To speed this
18	greg	2.4	* calculation, we sort the data into half-planes and apply
19			* simple tests to see which neighbor is closest in each
20			* direction. Once we have our approximate neighborhood
21	greg	2.3	* for a sample, we can use it in a modified Gaussian weighting
22	greg	2.4	* with allowing local anisotropy. Harmonic weighting is added
23	greg	2.3	* to reduce the influence of distant neighbors. This yields a
24			* smooth interpolation regardless of how the sample points are
25	greg	2.4	* initially distributed. Evaluation is accelerated by use of
26			* a fast approximation to the atan2(y,x) function.
27			****************************************************************/
28	greg	2.1
29			#include <stdio.h>
30			#include <stdlib.h>
31			#include "rtmath.h"
32			#include "interp2d.h"
33
34	greg	2.4	#define DECODE_DIA(ip,ed) ((ip)->dmin(1. + .5(ed)))
35			#define ENCODE_DIA(ip,d) ((int)(2.*(d)/(ip)->dmin) - 2)
36	greg	2.1
37			/* Sample order (private) */
38			typedef struct {
39			int si; /* sample index */
40			float dm; /* distance measure in this direction */
41			} SAMPORD;
42
43	greg	2.2	/* Allocate a new set of interpolation samples (caller assigns spt[] array) */
44	greg	2.1	INTERP2 *
45			interp2_alloc(int nsamps)
46			{
47			INTERP2 *nip;
48
49			if (nsamps <= 1)
50			return(NULL);
51
52			nip = (INTERP2 )malloc(sizeof(INTERP2) + sizeof(float)2*(nsamps-1));
53			if (nip == NULL)
54			return(NULL);
55
56			nip->ns = nsamps;
57	greg	2.4	nip->dmin = 1; /* default minimum diameter */
58	greg	2.1	nip->smf = NI2DSMF; /* default smoothing factor */
59	greg	2.4	nip->da = NULL;
60	greg	2.1	/* caller must assign spt[] array */
61			return(nip);
62			}
63
64	greg	2.2	/* Resize interpolation array (caller must assign any new values) */
65			INTERP2 *
66			interp2_realloc(INTERP2 *ip, int nsamps)
67			{
68			if (ip == NULL)
69			return(interp2_alloc(nsamps));
70			if (nsamps <= 1) {
71			interp2_free(ip);
72			return(NULL);
73			}
74			if (nsamps == ip->ns);
75			return(ip);
76	greg	2.4	if (ip->da != NULL) { /* will need to recompute distribution */
77			free(ip->da);
78			ip->da = NULL;
79	greg	2.2	}
80			ip = (INTERP2 )realloc(ip, sizeof(INTERP2)+sizeof(float)2*(nsamps-1));
81			if (ip == NULL)
82			return(NULL);
83			ip->ns = nsamps;
84			return(ip);
85			}
86
87	greg	2.5	/* Set minimum distance under which samples will start to merge */
88			void
89			interp2_spacing(INTERP2 *ip, double mind)
90			{
91			if (mind <= 0)
92			return;
93			if ((.998ip->dmin <= mind) && (mind <= 1.002ip->dmin))
94			return;
95			if (ip->da != NULL) { /* will need to recompute distribution */
96			free(ip->da);
97			ip->da = NULL;
98			}
99			ip->dmin = mind;
100			}
101
102			/* Modify smoothing parameter by the given factor */
103			void
104			interp2_smooth(INTERP2 *ip, double sf)
105			{
106			if ((ip->smf *= sf) < NI2DSMF)
107			ip->smf = NI2DSMF;
108			}
109
110	greg	2.1	/* private call-back to sort position index */
111			static int
112			cmp_spos(const void p1, const void p2)
113			{
114			const SAMPORD so1 = (const SAMPORD )p1;
115			const SAMPORD so2 = (const SAMPORD )p2;
116
117			if (so1->dm > so2->dm)
118			return 1;
119			if (so1->dm < so2->dm)
120			return -1;
121			return 0;
122			}
123
124	greg	2.2	/* private routine to order samples in a particular direction */
125			static void
126			sort_samples(SAMPORD sord, const INTERP2 ip, double ang)
127			{
128			const double cosd = cos(ang);
129			const double sind = sin(ang);
130			int i;
131
132			for (i = ip->ns; i--; ) {
133			sord[i].si = i;
134			sord[i].dm = cosdip->spt[i][0] + sindip->spt[i][1];
135			}
136			qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
137			}
138
139	greg	2.4	/* private routine to encode sample diameter with range checks */
140	greg	2.1	static int
141	greg	2.4	encode_diameter(const INTERP2 *ip, double d)
142	greg	2.1	{
143	greg	2.4	const int ed = ENCODE_DIA(ip, d);
144	greg	2.1
145	greg	2.4	if (ed <= 0)
146	greg	2.1	return(0);
147	greg	2.4	if (ed >= 0xffff)
148	greg	2.1	return(0xffff);
149	greg	2.4	return(ed);
150	greg	2.1	}
151
152	greg	2.2	/* (Re)compute anisotropic basis function interpolant (normally automatic) */
153			int
154			interp2_analyze(INTERP2 *ip)
155	greg	2.1	{
156			SAMPORD *sortord;
157	greg	2.2	int rightrndx, leftrndx, *endrndx;
158	greg	2.1	int bd;
159			/* sanity checks */
160	greg	2.4	if (ip == NULL \|\| (ip->ns <= 1) \| (ip->dmin <= 0))
161	greg	2.1	return(0);
162			/* need to allocate? */
163	greg	2.4	if (ip->da == NULL) {
164			ip->da = (unsigned short (*)[NI2DIR])malloc(
165	greg	2.1	sizeof(unsigned short)NI2DIRip->ns);
166	greg	2.4	if (ip->da == NULL)
167	greg	2.1	return(0);
168			}
169			/* get temporary arrays */
170			sortord = (SAMPORD )malloc(sizeof(SAMPORD)ip->ns);
171			rightrndx = (int )malloc(sizeof(int)ip->ns);
172			leftrndx = (int )malloc(sizeof(int)ip->ns);
173	greg	2.2	endrndx = (int )malloc(sizeof(int)ip->ns);
174			if ((sortord == NULL) \| (rightrndx == NULL) \|
175			(leftrndx == NULL) \| (endrndx == NULL))
176	greg	2.1	return(0);
177			/* run through bidirections */
178			for (bd = 0; bd < NI2DIR/2; bd++) {
179			const double ang = 2.PI/NI2DIRbd;
180	greg	2.2	int *sptr;
181	greg	2.1	int i;
182			/* create right reverse index */
183	greg	2.2	if (bd) { /* re-use from previous iteration? */
184			sptr = rightrndx;
185	greg	2.1	rightrndx = leftrndx;
186			leftrndx = sptr;
187	greg	2.2	} else { /* else sort first half-plane */
188			sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
189			for (i = ip->ns; i--; )
190	greg	2.1	rightrndx[sortord[i].si] = i;
191	greg	2.2	/* & store reverse order for later */
192			for (i = ip->ns; i--; )
193			endrndx[sortord[i].si] = ip->ns-1 - i;
194	greg	2.1	}
195			/* create new left reverse index */
196	greg	2.2	if (bd == NI2DIR/2 - 1) { /* use order from first iteration? */
197			sptr = leftrndx;
198			leftrndx = endrndx;
199			endrndx = sptr;
200			} else { /* else compute new half-plane */
201			sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
202			for (i = ip->ns; i--; )
203			leftrndx[sortord[i].si] = i;
204	greg	2.1	}
205			/* sort grid values in this direction */
206	greg	2.2	sort_samples(sortord, ip, ang);
207	greg	2.1	/* find nearest neighbors each side */
208	greg	2.2	for (i = ip->ns; i--; ) {
209	greg	2.3	const int ii = sortord[i].si;
210	greg	2.1	int j;
211	greg	2.3	/* preload with large radii */
212	greg	2.4	ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
213			.5*(sortord[ip->ns-1].dm - sortord[0].dm));
214	greg	2.1	for (j = i; ++j < ip->ns; ) /* nearest above */
215	greg	2.3	if (rightrndx[sortord[j].si] > rightrndx[ii] &&
216			leftrndx[sortord[j].si] < leftrndx[ii]) {
217	greg	2.4	ip->da[ii][bd] = encode_diameter(ip,
218			sortord[j].dm - sortord[i].dm);
219	greg	2.1	break;
220			}
221			for (j = i; j-- > 0; ) /* nearest below */
222	greg	2.3	if (rightrndx[sortord[j].si] < rightrndx[ii] &&
223			leftrndx[sortord[j].si] > leftrndx[ii]) {
224	greg	2.4	ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
225			sortord[i].dm - sortord[j].dm);
226	greg	2.1	break;
227			}
228			}
229			}
230			free(sortord); /* clean up */
231			free(rightrndx);
232			free(leftrndx);
233	greg	2.2	free(endrndx);
234	greg	2.1	return(1);
235			}
236
237	greg	2.3	/* private call returns raw weight for a particular sample */
238	greg	2.1	static double
239	greg	2.3	get_wt(const INTERP2 *ip, const int i, double x, double y)
240	greg	2.1	{
241	greg	2.3	double dir, rd, d2;
242	greg	2.1	int ri;
243			/* get relative direction */
244			x -= ip->spt[i][0];
245			y -= ip->spt[i][1];
246			dir = atan2a(y, x);
247			dir += 2.PI(dir < 0);
248			/* linear radius interpolation */
249			rd = dir * (NI2DIR/2./PI);
250			ri = (int)rd;
251			rd -= (double)ri;
252	greg	2.4	rd = (1.-rd)ip->da[i][ri] + rdip->da[i][(ri+1)%NI2DIR];
253			rd = ip->smf * DECODE_DIA(ip, rd);
254	greg	2.3	d2 = xx + yy;
255			/* Gaussian times harmonic weighting */
256	greg	2.4	return( exp(d2/(-2.rdrd)) * ip->dmin/(ip->dmin + sqrt(d2)) );
257	greg	2.1	}
258
259			/* Assign full set of normalized weights to interpolate the given position */
260			int
261			interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
262			{
263			double wnorm;
264			int i;
265
266			if ((wtv == NULL) \| (ip == NULL))
267			return(0);
268			/* need to compute interpolant? */
269	greg	2.4	if (ip->da == NULL && !interp2_analyze(ip))
270	greg	2.1	return(0);
271
272			wnorm = 0; /* compute raw weights */
273			for (i = ip->ns; i--; ) {
274	greg	2.3	double wt = get_wt(ip, i, x, y);
275	greg	2.1	wtv[i] = wt;
276			wnorm += wt;
277			}
278			if (wnorm <= 0) /* too far from all our samples! */
279			return(0);
280			wnorm = 1./wnorm; /* normalize weights */
281			for (i = ip->ns; i--; )
282			wtv[i] *= wnorm;
283			return(ip->ns); /* all done */
284			}
285
286
287			/* Get normalized weights and indexes for n best samples in descending order */
288			int
289			interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
290			{
291			int nn = 0;
292			double wnorm;
293			int i, j;
294
295			if ((n <= 0) \| (wt == NULL) \| (si == NULL) \| (ip == NULL))
296			return(0);
297			/* need to compute interpolant? */
298	greg	2.4	if (ip->da == NULL && !interp2_analyze(ip))
299	greg	2.1	return(0);
300			/* identify top n weights */
301			for (i = ip->ns; i--; ) {
302	greg	2.3	const double wti = get_wt(ip, i, x, y);
303	greg	2.1	for (j = nn; j > 0; j--) {
304	greg	2.3	if (wt[j-1] >= wti)
305	greg	2.1	break;
306			if (j < n) {
307			wt[j] = wt[j-1];
308			si[j] = si[j-1];
309			}
310			}
311			if (j < n) { /* add/insert sample */
312	greg	2.3	wt[j] = wti;
313	greg	2.1	si[j] = i;
314			nn += (nn < n);
315			}
316			}
317	greg	2.3	wnorm = 0; /* normalize sample weights */
318			for (j = nn; j--; )
319			wnorm += wt[j];
320	greg	2.1	if (wnorm <= 0)
321			return(0);
322			wnorm = 1./wnorm;
323			for (j = nn; j--; )
324			wt[j] *= wnorm;
325			return(nn); /* return actual # samples */
326			}
327
328			/* Free interpolant */
329			void
330			interp2_free(INTERP2 *ip)
331			{
332			if (ip == NULL)
333			return;
334	greg	2.4	if (ip->da != NULL)
335			free(ip->da);
336	greg	2.1	free(ip);
337			}