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.2 by greg, Sat Feb 9 17:39:21 2013 UTC vs.
Revision 2.13 by greg, Sat Feb 16 00:41:12 2013 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 + /* 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)
# 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->da = NULL;
74                                  /* caller must assign spt[] array */
75          return(nip);
76   }
# Line 70 | Line 85 | interp2_realloc(INTERP2 *ip, int nsamps)
85                  interp2_free(ip);
86                  return(NULL);
87          }
88 <        if (nsamps == ip->ns);
88 >        if (nsamps == ip->ns)
89                  return(ip);
90 <        if (ip->ra != NULL) {   /* will need to recompute distribution */
91 <                free(ip->ra);
92 <                ip->ra = NULL;
90 >        if (ip->da != NULL) {   /* will need to recompute distribution */
91 >                free(ip->da);
92 >                ip->da = NULL;
93          }
94          ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
95          if (ip == NULL)
# Line 83 | Line 98 | interp2_realloc(INTERP2 *ip, int nsamps)
98          return(ip);
99   }
100  
101 + /* Set minimum distance under which samples will start to merge */
102 + void
103 + interp2_spacing(INTERP2 *ip, double mind)
104 + {
105 +        if (mind <= 0)
106 +                return;
107 +        if ((.998*ip->dmin <= mind) & (mind <= 1.002*ip->dmin))
108 +                return;
109 +        if (ip->da != NULL) {   /* will need to recompute distribution */
110 +                free(ip->da);
111 +                ip->da = NULL;
112 +        }
113 +        ip->dmin = mind;
114 + }
115 +
116 + /* Compute unnormalized weight for a position relative to a sample */
117 + double
118 + interp2_wti(INTERP2 *ip, const int i, double x, double y)
119 + {
120 +        double          dir, rd, r2, d2;
121 +        int             ri;
122 +                                /* get relative direction */
123 +        x -= ip->spt[i][0];
124 +        y -= ip->spt[i][1];
125 +        dir = atan2a(y, x);
126 +        dir += 2.*PI*(dir < 0);
127 +                                /* linear radius interpolation */
128 +        rd = dir * (NI2DIR/2./PI);
129 +        ri = (int)rd;
130 +        rd -= (double)ri;
131 +        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
132 +        rd = ip->smf * DECODE_DIA(ip, rd);
133 +        r2 = 2.*rd*rd;
134 +        d2 = x*x + y*y;
135 +        if (d2 > 21.*r2)        /* result would be < 1e-9 */
136 +                return(.0);
137 +                                /* Gaussian times harmonic weighting */
138 +        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
139 + }
140 +
141 + /* private call to get grid flag index */
142 + static int
143 + interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
144 + {
145 +        int     inbounds = 0;
146 +
147 +        if (ip == NULL)         /* paranoia */
148 +                return(-1);
149 +                                /* need to compute interpolant? */
150 +        if (ip->da == NULL && !interp2_analyze(ip))
151 +                return(-1);
152 +                                /* get x & y grid positions */
153 +        fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
154 +
155 +        if (fgi[0] >= NI2DIM)
156 +                fgi[0] = NI2DIM-1;
157 +        else if (fgi[0] < 0)
158 +                fgi[0] = 0;
159 +        else
160 +                ++inbounds;
161 +
162 +        fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
163 +
164 +        if (fgi[1] >= NI2DIM)
165 +                fgi[1] = NI2DIM-1;
166 +        else if (fgi[1] < 0)
167 +                fgi[1] = 0;
168 +        else
169 +                ++inbounds;
170 +
171 +        return(inbounds == 2);
172 + }
173 +
174 + #define setflg(fl,xi,yi)        ((fl)[yi] |= 1<<(xi))
175 +
176 + #define chkflg(fl,xi,yi)        ((fl)[yi]>>(xi) & 1)
177 +
178 + /* private flood function to determine sample influence */
179 + static void
180 + influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
181 +                        int xfi, int yfi)
182 + {
183 +        double  gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
184 +                                        ip->smin[0];
185 +        double  gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
186 +                                        ip->smin[1];
187 +        double  dx = gx - ip->spt[i][0];
188 +        double  dy = gy - ip->spt[i][1];
189 +
190 +        setflg(visited, xfi, yfi);
191 +
192 +        if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
193 +                return;
194 +
195 +        setflg(ip->da[i].infl, xfi, yfi);
196 +
197 +        if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
198 +                influence_flood(ip, i, visited, xfi-1, yfi);
199 +
200 +        if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
201 +                influence_flood(ip, i, visited, xfi, yfi-1);
202 +
203 +        if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
204 +                influence_flood(ip, i, visited, xfi+1, yfi);
205 +
206 +        if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
207 +                influence_flood(ip, i, visited, xfi, yfi+1);
208 + }
209 +
210 + /* private call to compute sample influence maps */
211 + static void
212 + map_influence(INTERP2 *ip)
213 + {
214 +        unsigned short  visited[NI2DIM];
215 +        int             fgi[2];
216 +        int             i, j;
217 +
218 +        for (i = ip->ns; i--; ) {
219 +                for (j = NI2DIM; j--; ) {
220 +                        ip->da[i].infl[j] = 0;
221 +                        visited[j] = 0;
222 +                }
223 +                interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
224 +
225 +                influence_flood(ip, i, visited, fgi[0], fgi[1]);
226 +        }
227 + }
228 +
229 + /* Modify smoothing parameter by the given factor */
230 + void
231 + interp2_smooth(INTERP2 *ip, double sf)
232 + {
233 +        float   old_smf = ip->smf;
234 +
235 +        if ((ip->smf *= sf) < NI2DSMF)
236 +                ip->smf = NI2DSMF;
237 +                                        /* need to recompute influence maps? */
238 +        if (ip->da != NULL && (old_smf*.85 > ip->smf) |
239 +                                (ip->smf > old_smf*1.15))
240 +                map_influence(ip);
241 + }
242 +
243   /* private call-back to sort position index */
244   static int
245   cmp_spos(const void *p1, const void *p2)
# Line 112 | Line 269 | sort_samples(SAMPORD *sord, const INTERP2 *ip, double
269          qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
270   }
271  
115 – /* private routine to encode radius with range checks */
116 – static int
117 – encode_radius(const INTERP2 *ip, double r)
118 – {
119 –        const int       er = ENCODE_RAD(ip, r);
120 –
121 –        if (er <= 0)
122 –                return(0);
123 –        if (er >= 0xffff)
124 –                return(0xffff);
125 –        return(er);
126 – }
127 –
272   /* (Re)compute anisotropic basis function interpolant (normally automatic) */
273   int
274   interp2_analyze(INTERP2 *ip)
275   {
276          SAMPORD *sortord;
277          int     *rightrndx, *leftrndx, *endrndx;
278 <        int     bd;
278 >        int     i, bd;
279                                          /* sanity checks */
280 <        if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
280 >        if (ip == NULL)
281                  return(0);
282 <                                        /* need to allocate? */
283 <        if (ip->ra == NULL) {
284 <                ip->ra = (unsigned short (*)[NI2DIR])malloc(
141 <                                sizeof(unsigned short)*NI2DIR*ip->ns);
142 <                if (ip->ra == NULL)
143 <                        return(0);
282 >        if (ip->da != NULL) {           /* free previous data if any */
283 >                free(ip->da);
284 >                ip->da = NULL;
285          }
286 <                                        /* get temporary arrays */
286 >        if ((ip->ns <= 1) | (ip->dmin <= 0))
287 >                return(0);
288 >                                        /* compute sample domain */
289 >        ip->smin[0] = ip->smin[1] = FHUGE;
290 >        ip->smax[0] = ip->smax[1] = -FHUGE;
291 >        for (i = ip->ns; i--; ) {
292 >                if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
293 >                if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
294 >                if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
295 >                if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
296 >        }
297 >        ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
298 >                        (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
299 >                                (1./NI2DIM/NI2DIM);
300 >        if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
301 >                return(0);
302 >                                        /* allocate analysis data */
303 >        ip->da = (struct interp2_samp *)malloc(
304 >                                sizeof(struct interp2_samp)*ip->ns );
305 >        if (ip->da == NULL)
306 >                return(0);
307 >                                        /* allocate temporary arrays */
308          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
309          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
310          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
# Line 154 | Line 316 | interp2_analyze(INTERP2 *ip)
316          for (bd = 0; bd < NI2DIR/2; bd++) {
317              const double        ang = 2.*PI/NI2DIR*bd;
318              int                 *sptr;
157 –            int                 i;
319                                          /* create right reverse index */
320              if (bd) {                   /* re-use from previous iteration? */
321                  sptr = rightrndx;
# Line 182 | Line 343 | interp2_analyze(INTERP2 *ip)
343              sort_samples(sortord, ip, ang);
344                                          /* find nearest neighbors each side */
345              for (i = ip->ns; i--; ) {
346 <                const int       rpos = rightrndx[sortord[i].si];
186 <                const int       lpos = leftrndx[sortord[i].si];
346 >                const int       ii = sortord[i].si;
347                  int             j;
348 <                                        /* preload with large radius */
349 <                ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
350 <                            .25*(sortord[ip->ns-1].dm - sortord[0].dm));
348 >                                        /* preload with large radii */
349 >                ip->da[ii].dia[bd] =
350 >                ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
351 >                                .5*(sortord[ip->ns-1].dm - sortord[0].dm));
352                  for (j = i; ++j < ip->ns; )     /* nearest above */
353 <                    if (rightrndx[sortord[j].si] > rpos &&
354 <                                    leftrndx[sortord[j].si] < lpos) {
355 <                        ip->ra[i][bd] = encode_radius(ip,
356 <                                        .5*(sortord[j].dm - sortord[i].dm));
353 >                    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
354 >                                    leftrndx[sortord[j].si] < leftrndx[ii]) {
355 >                        ip->da[ii].dia[bd] = encode_diameter(ip,
356 >                                                sortord[j].dm - sortord[i].dm);
357                          break;
358                      }
359                  for (j = i; j-- > 0; )          /* nearest below */
360 <                    if (rightrndx[sortord[j].si] < rpos &&
361 <                                    leftrndx[sortord[j].si] > lpos) {
362 <                        ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
363 <                                        .5*(sortord[i].dm - sortord[j].dm));
360 >                    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
361 >                                    leftrndx[sortord[j].si] > leftrndx[ii]) {
362 >                        ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
363 >                                                sortord[i].dm - sortord[j].dm);
364                          break;
365                      }
366              }
367          }
368 <        free(sortord);                  /* clean up */
368 >        free(sortord);                  /* release temp arrays */
369          free(rightrndx);
370          free(leftrndx);
371          free(endrndx);
372 <        return(1);
372 >                                        /* map sample influence areas */
373 >        map_influence(ip);
374 >        return(1);                      /* all done */
375   }
376  
214 – /* private call returns log of raw weight for a particular sample */
215 – static double
216 – get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
217 – {
218 –        double  dir, rd;
219 –        int     ri;
220 –                                /* get relative direction */
221 –        x -= ip->spt[i][0];
222 –        y -= ip->spt[i][1];
223 –        dir = atan2a(y, x);
224 –        dir += 2.*PI*(dir < 0);
225 –                                /* linear radius interpolation */
226 –        rd = dir * (NI2DIR/2./PI);
227 –        ri = (int)rd;
228 –        rd -= (double)ri;
229 –        rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
230 –        rd = ip->smf * DECODE_RAD(ip, rd);
231 –                                /* return log of Gaussian weight */
232 –        return( (x*x + y*y) / (-2.*rd*rd) );
233 – }
234 –
377   /* Assign full set of normalized weights to interpolate the given position */
378   int
379   interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
380   {
381          double  wnorm;
382 +        int     fgi[2];
383          int     i;
384  
385 <        if ((wtv == NULL) | (ip == NULL))
385 >        if (wtv == NULL)
386                  return(0);
387 <                                        /* need to compute interpolant? */
388 <        if (ip->ra == NULL && !interp2_analyze(ip))
387 >                                        /* get flag position */
388 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
389                  return(0);
390  
391          wnorm = 0;                      /* compute raw weights */
392 <        for (i = ip->ns; i--; ) {
393 <                double  wt = get_ln_wt(ip, i, x, y);
394 <                if (wt < -21.) {
252 <                        wtv[i] = 0;     /* ignore weights < 1e-9 */
253 <                        continue;
254 <                }
255 <                wt = exp(wt);           /* Gaussian weight */
392 >        for (i = ip->ns; i--; )
393 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
394 >                double  wt = interp2_wti(ip, i, x, y);
395                  wtv[i] = wt;
396                  wnorm += wt;
397 <        }
397 >            } else
398 >                wtv[i] = 0;
399          if (wnorm <= 0)                 /* too far from all our samples! */
400                  return(0);
401          wnorm = 1./wnorm;               /* normalize weights */
# Line 270 | Line 410 | int
410   interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
411   {
412          int     nn = 0;
413 +        int     fgi[2];
414          double  wnorm;
415          int     i, j;
416  
417 <        if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
417 >        if ((n <= 0) | (wt == NULL) | (si == NULL))
418                  return(0);
419 <                                        /* need to compute interpolant? */
420 <        if (ip->ra == NULL && !interp2_analyze(ip))
419 >                                        /* get flag position */
420 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
421                  return(0);
422                                          /* identify top n weights */
423 <        for (i = ip->ns; i--; ) {
424 <                const double    lnwt = get_ln_wt(ip, i, x, y);
423 >        for (i = ip->ns; i--; )
424 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
425 >                const double    wti = interp2_wti(ip, i, x, y);
426                  for (j = nn; j > 0; j--) {
427 <                        if (wt[j-1] >= lnwt)
427 >                        if (wt[j-1] >= wti)
428                                  break;
429                          if (j < n) {
430                                  wt[j] = wt[j-1];
# Line 290 | Line 432 | interp2_topsamp(float wt[], int si[], const int n, INT
432                          }
433                  }
434                  if (j < n) {            /* add/insert sample */
435 <                        wt[j] = lnwt;
435 >                        wt[j] = wti;
436                          si[j] = i;
437                          nn += (nn < n);
438                  }
439 <        }
440 <        wnorm = 0;                      /* exponentiate and normalize */
441 <        for (j = nn; j--; ) {
442 <                double  dwt = exp(wt[j]);
301 <                wt[j] = dwt;
302 <                wnorm += dwt;
303 <        }
439 >            }
440 >        wnorm = 0;                      /* normalize sample weights */
441 >        for (j = nn; j--; )
442 >                wnorm += wt[j];
443          if (wnorm <= 0)
444                  return(0);
445          wnorm = 1./wnorm;
# Line 315 | Line 454 | interp2_free(INTERP2 *ip)
454   {
455          if (ip == NULL)
456                  return;
457 <        if (ip->ra != NULL)
458 <                free(ip->ra);
457 >        if (ip->da != NULL)
458 >                free(ip->da);
459          free(ip);
460   }

Diff Legend

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