ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/interp2d.c
(Generate patch)

Comparing ray/src/common/interp2d.c (file contents):
Revision 2.1 by greg, Sat Feb 9 00:55:40 2013 UTC vs.
Revision 2.17 by greg, Fri Mar 19 21:16:15 2021 UTC

# Line 9 | Line 9 | static const char RCSid[] = "$Id$";
9  
10   #include "copyright.h"
11  
12 < /*************************************************************
12 > /***************************************************************
13   * 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 < * calculation, we sort the data into 3 half-planes and
19 < * perform simple tests to see which neighbor is closest in
20 < * a each direction.  Once we have our approximate neighborhood
21 < * for a sample, we can use it in a Gaussian weighting scheme
22 < * with anisotropic surround.  This gives us a fairly smooth
23 < * interpolation however the sample points may be initially
24 < * distributed.  Evaluation is accelerated by use of a fast
25 < * approximation to the atan2(y,x) function.
26 < **************************************************************/
18 > * calculation, we sort the data into half-planes and apply
19 > * simple tests to see which neighbor is closest in each
20 > * angular slice.  Once we have our approximate neighborhood
21 > * for a sample, we can use it in a modified Gaussian weighting
22 > * with allowing local anisotropy.  Harmonic weighting is added
23 > * to reduce the influence of distant neighbors.  This yields a
24 > * smooth interpolation regardless of how the sample points are
25 > * initially distributed.  Evaluation is accelerated by use of
26 > * a fast approximation to the atan2(y,x) function and a low-res
27 > * map indicating where sample weights are significant.
28 > ****************************************************************/
29  
30   #include <stdio.h>
31   #include <stdlib.h>
32   #include "rtmath.h"
33   #include "interp2d.h"
34  
35 < #define DECODE_RAD(ip,er)       ((ip)->rmin*(1. + .5*(er)))
36 < #define ENCODE_RAD(ip,r)        ((int)(2.*(r)/(ip)->rmin) - 2)
35 > #define DECODE_DIA(ip,ed)       ((ip)->dmin*(1. + .5*(ed)))
36 > #define ENCODE_DIA(ip,d)        ((int)(2.*(d)/(ip)->dmin) - 2)
37  
38   /* Sample order (private) */
39   typedef struct {
# Line 39 | Line 41 | typedef struct {
41          float   dm;             /* distance measure in this direction */
42   } SAMPORD;
43  
44 < /* Allocate a new set of interpolation samples */
44 > /* private routine to encode sample diameter with range checks */
45 > static int
46 > encode_diameter(const INTERP2 *ip, double d)
47 > {
48 >        const int       ed = ENCODE_DIA(ip, d);
49 >
50 >        if (ed <= 0)
51 >                return(0);
52 >        if (ed >= 0xffff)
53 >                return(0xffff);
54 >        return(ed);
55 > }
56 >
57 > /* Allocate a new set of interpolation samples (caller assigns spt[] array) */
58   INTERP2 *
59   interp2_alloc(int nsamps)
60   {
# Line 53 | Line 68 | interp2_alloc(int nsamps)
68                  return(NULL);
69  
70          nip->ns = nsamps;
71 <        nip->rmin = .5;         /* default radius minimum */
71 >        nip->dmin = 1;          /* default minimum diameter */
72          nip->smf = NI2DSMF;     /* default smoothing factor */
73 <        nip->ra = NULL;
73 >        nip->c_data = NULL;
74 >        nip->da = NULL;
75                                  /* caller must assign spt[] array */
76          return(nip);
77   }
78  
79 + /* Resize interpolation array (caller must assign any new values) */
80 + INTERP2 *
81 + interp2_realloc(INTERP2 *ip, int nsamps)
82 + {
83 +        INTERP2 *old_ip = ip;
84 +
85 +        if (ip == NULL)
86 +                return(interp2_alloc(nsamps));
87 +        if (nsamps <= 1) {
88 +                interp2_free(ip);
89 +                return(NULL);
90 +        }
91 +        if (nsamps == ip->ns)
92 +                return(ip);
93 +        if (ip->da != NULL) {   /* will need to recompute distribution */
94 +                free(ip->da);
95 +                ip->da = NULL;
96 +        }
97 +        ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
98 +        if (ip == NULL) {
99 +                if (nsamps <= ip->ns) {
100 +                        ip = old_ip;
101 +                } else {
102 +                        free(old_ip);
103 +                        return(NULL);
104 +                }
105 +        }
106 +        ip->ns = nsamps;
107 +        return(ip);
108 + }
109 +
110 + /* Set minimum distance under which samples will start to merge */
111 + void
112 + interp2_spacing(INTERP2 *ip, double mind)
113 + {
114 +        if (mind <= 0)
115 +                return;
116 +        if ((.998*ip->dmin <= mind) & (mind <= 1.002*ip->dmin))
117 +                return;
118 +        if (ip->da != NULL) {   /* will need to recompute distribution */
119 +                free(ip->da);
120 +                ip->da = NULL;
121 +        }
122 +        ip->dmin = mind;
123 + }
124 +
125 + /* Compute unnormalized weight for a position relative to a sample */
126 + double
127 + interp2_wti(INTERP2 *ip, const int i, double x, double y)
128 + {
129 +        double          dir, rd, r2, d2;
130 +        int             ri;
131 +                                /* get relative direction */
132 +        x -= ip->spt[i][0];
133 +        y -= ip->spt[i][1];
134 +        dir = atan2a(y, x);
135 +        dir += 2.*PI*(dir < 0);
136 +                                /* linear radius interpolation */
137 +        rd = dir * (NI2DIR/2./PI);
138 +        ri = (int)rd;
139 +        rd -= (double)ri;
140 +        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
141 +        rd = ip->smf * DECODE_DIA(ip, rd);
142 +        r2 = 2.*rd*rd;
143 +        d2 = x*x + y*y;
144 +        if (d2 > 21.*r2)        /* result would be < 1e-9 */
145 +                return(.0);
146 +                                /* Gaussian times harmonic weighting */
147 +        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
148 + }
149 +
150 + /* private call to get grid flag index */
151 + static int
152 + interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
153 + {
154 +        int     inbounds = 0;
155 +
156 +        if (ip == NULL)         /* paranoia */
157 +                return(-1);
158 +                                /* need to compute interpolant? */
159 +        if (ip->da == NULL && !interp2_analyze(ip))
160 +                return(-1);
161 +                                /* get x & y grid positions */
162 +        fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
163 +
164 +        if (fgi[0] >= NI2DIM)
165 +                fgi[0] = NI2DIM-1;
166 +        else if (fgi[0] < 0)
167 +                fgi[0] = 0;
168 +        else
169 +                ++inbounds;
170 +
171 +        fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
172 +
173 +        if (fgi[1] >= NI2DIM)
174 +                fgi[1] = NI2DIM-1;
175 +        else if (fgi[1] < 0)
176 +                fgi[1] = 0;
177 +        else
178 +                ++inbounds;
179 +
180 +        return(inbounds == 2);
181 + }
182 +
183 + #define setflg(fl,xi,yi)        ((fl)[yi] |= 1<<(xi))
184 +
185 + #define chkflg(fl,xi,yi)        ((fl)[yi]>>(xi) & 1)
186 +
187 + /* private flood function to determine sample influence */
188 + static void
189 + influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
190 +                        int xfi, int yfi)
191 + {
192 +        double  gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
193 +                                        ip->smin[0];
194 +        double  gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
195 +                                        ip->smin[1];
196 +        double  dx = gx - ip->spt[i][0];
197 +        double  dy = gy - ip->spt[i][1];
198 +
199 +        setflg(visited, xfi, yfi);
200 +
201 +        if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
202 +                return;
203 +
204 +        setflg(ip->da[i].infl, xfi, yfi);
205 +
206 +        if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
207 +                influence_flood(ip, i, visited, xfi-1, yfi);
208 +
209 +        if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
210 +                influence_flood(ip, i, visited, xfi, yfi-1);
211 +
212 +        if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
213 +                influence_flood(ip, i, visited, xfi+1, yfi);
214 +
215 +        if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
216 +                influence_flood(ip, i, visited, xfi, yfi+1);
217 + }
218 +
219 + /* private call to compute sample influence maps */
220 + static void
221 + map_influence(INTERP2 *ip)
222 + {
223 +        unsigned short  visited[NI2DIM];
224 +        int             fgi[2];
225 +        int             i, j;
226 +
227 +        for (i = ip->ns; i--; ) {
228 +                for (j = NI2DIM; j--; ) {
229 +                        ip->da[i].infl[j] = 0;
230 +                        visited[j] = 0;
231 +                }
232 +                interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
233 +
234 +                influence_flood(ip, i, visited, fgi[0], fgi[1]);
235 +        }
236 + }
237 +
238 + /* Modify smoothing parameter by the given factor */
239 + void
240 + interp2_smooth(INTERP2 *ip, double sf)
241 + {
242 +        float   old_smf = ip->smf;
243 +
244 +        if ((ip->smf *= sf) < NI2DSMF)
245 +                ip->smf = NI2DSMF;
246 +                                        /* need to recompute influence maps? */
247 +        if (ip->da != NULL && (old_smf*.85 > ip->smf) |
248 +                                (ip->smf > old_smf*1.15))
249 +                map_influence(ip);
250 + }
251 +
252   /* private call-back to sort position index */
253   static int
254   cmp_spos(const void *p1, const void *p2)
# Line 74 | Line 263 | cmp_spos(const void *p1, const void *p2)
263          return 0;
264   }
265  
266 < /* private routine to encode radius with range checks */
267 < static int
268 < encode_radius(const INTERP2 *ip, double r)
266 > /* private routine to order samples in a particular direction */
267 > static void
268 > sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang)
269   {
270 <        const int       er = ENCODE_RAD(ip, r);
270 >        const double    cosd = cos(ang);
271 >        const double    sind = sin(ang);
272 >        int             i;
273  
274 <        if (er <= 0)
275 <                return(0);
276 <        if (er >= 0xffff)
277 <                return(0xffff);
278 <        return(er);
274 >        for (i = ip->ns; i--; ) {
275 >                sord[i].si = i;
276 >                sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
277 >        }
278 >        qsort(sord, ip->ns, sizeof(SAMPORD), cmp_spos);
279   }
280  
281 < /* Compute anisotropic Gaussian basis function interpolant */
282 < static int
283 < interp2_compute(INTERP2 *ip)
281 > /* (Re)compute anisotropic basis function interpolant (normally automatic) */
282 > int
283 > interp2_analyze(INTERP2 *ip)
284   {
285          SAMPORD *sortord;
286 <        int     *rightrndx, *leftrndx;
287 <        int     bd;
286 >        int     *rightrndx, *leftrndx, *endrndx;
287 >        int     i, bd;
288                                          /* sanity checks */
289 <        if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
289 >        if (ip == NULL)
290                  return(0);
291 <                                        /* need to allocate? */
292 <        if (ip->ra == NULL) {
293 <                ip->ra = (unsigned short (*)[NI2DIR])malloc(
103 <                                sizeof(unsigned short)*NI2DIR*ip->ns);
104 <                if (ip->ra == NULL)
105 <                        return(0);
291 >        if (ip->da != NULL) {           /* free previous data if any */
292 >                free(ip->da);
293 >                ip->da = NULL;
294          }
295 <                                        /* get temporary arrays */
295 >        if ((ip->ns <= 1) | (ip->dmin <= 0))
296 >                return(0);
297 >                                        /* compute sample domain */
298 >        ip->smin[0] = ip->smin[1] = FHUGE;
299 >        ip->smax[0] = ip->smax[1] = -FHUGE;
300 >        for (i = ip->ns; i--; ) {
301 >                if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
302 >                if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
303 >                if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
304 >                if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
305 >        }
306 >        ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
307 >                        (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
308 >                                (1./NI2DIM/NI2DIM);
309 >        if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
310 >                return(0);
311 >                                        /* allocate analysis data */
312 >        ip->da = (struct interp2_samp *)malloc(
313 >                                sizeof(struct interp2_samp)*ip->ns );
314 >        if (ip->da == NULL)
315 >                return(0);
316 >                                        /* allocate temporary arrays */
317          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
318          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
319          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
320 <        if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL))
320 >        endrndx = (int *)malloc(sizeof(int)*ip->ns);
321 >        if ((sortord == NULL) | (rightrndx == NULL) |
322 >                        (leftrndx == NULL) | (endrndx == NULL))
323                  return(0);
324                                          /* run through bidirections */
325          for (bd = 0; bd < NI2DIR/2; bd++) {
326              const double        ang = 2.*PI/NI2DIR*bd;
327 <            double              cosd, sind;
117 <            int                 i;
327 >            int                 *sptr;
328                                          /* create right reverse index */
329 <            if (bd) {                   /* re-use from prev. iteration? */
330 <                int     *sptr = rightrndx;
329 >            if (bd) {                   /* re-use from previous iteration? */
330 >                sptr = rightrndx;
331                  rightrndx = leftrndx;
332                  leftrndx = sptr;
333 <            } else {                    /* else compute it */
334 <                cosd = cos(ang + (PI/2. - PI/NI2DIR));
335 <                sind = sin(ang + (PI/2. - PI/NI2DIR));
126 <                for (i = 0; i < ip->ns; i++) {
127 <                    sortord[i].si = i;
128 <                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
129 <                }
130 <                qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
131 <                for (i = 0; i < ip->ns; i++)
333 >            } else {                    /* else sort first half-plane */
334 >                sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
335 >                for (i = ip->ns; i--; )
336                      rightrndx[sortord[i].si] = i;
337 +                                        /* & store reverse order for later */
338 +                for (i = ip->ns; i--; )
339 +                    endrndx[sortord[i].si] = ip->ns-1 - i;
340              }
341                                          /* create new left reverse index */
342 <            cosd = cos(ang + (PI/2. + PI/NI2DIR));
343 <            sind = sin(ang + (PI/2. + PI/NI2DIR));
344 <            for (i = 0; i < ip->ns; i++) {
345 <                    sortord[i].si = i;
346 <                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
347 <            }
348 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
142 <            for (i = 0; i < ip->ns; i++)
342 >            if (bd == NI2DIR/2 - 1) {   /* use order from first iteration? */
343 >                sptr = leftrndx;
344 >                leftrndx = endrndx;
345 >                endrndx = sptr;
346 >            } else {                    /* else compute new half-plane */
347 >                sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
348 >                for (i = ip->ns; i--; )
349                      leftrndx[sortord[i].si] = i;
144 –                                        /* sort grid values in this direction */
145 –            cosd = cos(ang);
146 –            sind = sin(ang);
147 –            for (i = 0; i < ip->ns; i++) {
148 –                    sortord[i].si = i;
149 –                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
350              }
351 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
351 >                                        /* sort grid values in this direction */
352 >            sort_samples(sortord, ip, ang);
353                                          /* find nearest neighbors each side */
354 <            for (i = 0; i < ip->ns; i++) {
355 <                const int       rpos = rightrndx[sortord[i].si];
155 <                const int       lpos = leftrndx[sortord[i].si];
354 >            for (i = ip->ns; i--; ) {
355 >                const int       ii = sortord[i].si;
356                  int             j;
357 <                                        /* preload with large radius */
358 <                ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
359 <                            .25*(sortord[ip->ns-1].dm - sortord[0].dm));
357 >                                        /* preload with large radii */
358 >                ip->da[ii].dia[bd] =
359 >                ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
360 >                                .5*(sortord[ip->ns-1].dm - sortord[0].dm));
361                  for (j = i; ++j < ip->ns; )     /* nearest above */
362 <                    if (rightrndx[sortord[j].si] > rpos &&
363 <                                    leftrndx[sortord[j].si] < lpos) {
364 <                        ip->ra[i][bd] = encode_radius(ip,
365 <                                        .5*(sortord[j].dm - sortord[i].dm));
362 >                    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
363 >                                    leftrndx[sortord[j].si] < leftrndx[ii]) {
364 >                        ip->da[ii].dia[bd] = encode_diameter(ip,
365 >                                                sortord[j].dm - sortord[i].dm);
366                          break;
367                      }
368                  for (j = i; j-- > 0; )          /* nearest below */
369 <                    if (rightrndx[sortord[j].si] < rpos &&
370 <                                    leftrndx[sortord[j].si] > lpos) {
371 <                        ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
372 <                                        .5*(sortord[i].dm - sortord[j].dm));
369 >                    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
370 >                                    leftrndx[sortord[j].si] > leftrndx[ii]) {
371 >                        ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
372 >                                                sortord[i].dm - sortord[j].dm);
373                          break;
374                      }
375              }
376          }
377 <        free(sortord);                  /* clean up */
377 >        free(sortord);                  /* release temp arrays */
378          free(rightrndx);
379          free(leftrndx);
380 <        return(1);
380 >        free(endrndx);
381 >                                        /* map sample influence areas */
382 >        map_influence(ip);
383 >        return(1);                      /* all done */
384   }
385  
182 – /* private call returns log of raw weight for a particular sample */
183 – static double
184 – get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
185 – {
186 –        double  dir, rd;
187 –        int     ri;
188 –                                /* get relative direction */
189 –        x -= ip->spt[i][0];
190 –        y -= ip->spt[i][1];
191 –        dir = atan2a(y, x);
192 –        dir += 2.*PI*(dir < 0);
193 –                                /* linear radius interpolation */
194 –        rd = dir * (NI2DIR/2./PI);
195 –        ri = (int)rd;
196 –        rd -= (double)ri;
197 –        rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
198 –        rd = ip->smf * DECODE_RAD(ip, rd);
199 –                                /* return log of Gaussian weight */
200 –        return( (x*x + y*y) / (-2.*rd*rd) );
201 – }
202 –
386   /* Assign full set of normalized weights to interpolate the given position */
387   int
388   interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
389   {
390          double  wnorm;
391 +        int     fgi[2];
392          int     i;
393  
394 <        if ((wtv == NULL) | (ip == NULL))
394 >        if (wtv == NULL)
395                  return(0);
396 <                                        /* need to compute interpolant? */
397 <        if (ip->ra == NULL && !interp2_compute(ip))
396 >                                        /* get flag position */
397 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
398                  return(0);
399  
400          wnorm = 0;                      /* compute raw weights */
401 <        for (i = ip->ns; i--; ) {
402 <                double  wt = get_ln_wt(ip, i, x, y);
403 <                if (wt < -21.) {
220 <                        wtv[i] = 0;     /* ignore weights < 1e-9 */
221 <                        continue;
222 <                }
223 <                wt = exp(wt);           /* Gaussian weight */
401 >        for (i = ip->ns; i--; )
402 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
403 >                double  wt = interp2_wti(ip, i, x, y);
404                  wtv[i] = wt;
405                  wnorm += wt;
406 <        }
406 >            } else
407 >                wtv[i] = 0;
408          if (wnorm <= 0)                 /* too far from all our samples! */
409                  return(0);
410          wnorm = 1./wnorm;               /* normalize weights */
# Line 238 | Line 419 | int
419   interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
420   {
421          int     nn = 0;
422 +        int     fgi[2];
423          double  wnorm;
424          int     i, j;
425  
426 <        if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
426 >        if ((n <= 0) | (wt == NULL) | (si == NULL))
427                  return(0);
428 <                                        /* need to compute interpolant? */
429 <        if (ip->ra == NULL && !interp2_compute(ip))
428 >                                        /* get flag position */
429 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
430                  return(0);
431                                          /* identify top n weights */
432 <        for (i = ip->ns; i--; ) {
433 <                const double    lnwt = get_ln_wt(ip, i, x, y);
432 >        for (i = ip->ns; i--; )
433 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
434 >                const double    wti = interp2_wti(ip, i, x, y);
435                  for (j = nn; j > 0; j--) {
436 <                        if (wt[j-1] >= lnwt)
436 >                        if (wt[j-1] >= wti)
437                                  break;
438                          if (j < n) {
439                                  wt[j] = wt[j-1];
# Line 258 | Line 441 | interp2_topsamp(float wt[], int si[], const int n, INT
441                          }
442                  }
443                  if (j < n) {            /* add/insert sample */
444 <                        wt[j] = lnwt;
444 >                        wt[j] = wti;
445                          si[j] = i;
446                          nn += (nn < n);
447                  }
448 <        }
449 <        wnorm = 0;                      /* exponentiate and normalize */
450 <        for (j = nn; j--; ) {
451 <                double  dwt = exp(wt[j]);
269 <                wt[j] = dwt;
270 <                wnorm += dwt;
271 <        }
448 >            }
449 >        wnorm = 0;                      /* normalize sample weights */
450 >        for (j = nn; j--; )
451 >                wnorm += wt[j];
452          if (wnorm <= 0)
453                  return(0);
454          wnorm = 1./wnorm;
# Line 283 | Line 463 | interp2_free(INTERP2 *ip)
463   {
464          if (ip == NULL)
465                  return;
466 <        if (ip->ra != NULL)
467 <                free(ip->ra);
466 >        if (ip->da != NULL)
467 >                free(ip->da);
468          free(ip);
469   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines