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.11 by greg, Fri Feb 15 01:26:47 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 an array
27 > * of flags indicating where weights are (nearly) zero.
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 > /* Allocate a new set of interpolation samples (caller assigns spt[] array) */
45   INTERP2 *
46   interp2_alloc(int nsamps)
47   {
# Line 53 | Line 55 | interp2_alloc(int nsamps)
55                  return(NULL);
56  
57          nip->ns = nsamps;
58 <        nip->rmin = .5;         /* default radius minimum */
58 >        nip->dmin = 1;          /* default minimum diameter */
59          nip->smf = NI2DSMF;     /* default smoothing factor */
60 <        nip->ra = NULL;
60 >        nip->da = NULL;
61                                  /* caller must assign spt[] array */
62          return(nip);
63   }
64  
65 + /* Resize interpolation array (caller must assign any new values) */
66 + INTERP2 *
67 + interp2_realloc(INTERP2 *ip, int nsamps)
68 + {
69 +        if (ip == NULL)
70 +                return(interp2_alloc(nsamps));
71 +        if (nsamps <= 1) {
72 +                interp2_free(ip);
73 +                return(NULL);
74 +        }
75 +        if (nsamps == ip->ns)
76 +                return(ip);
77 +        if (ip->da != NULL) {   /* will need to recompute distribution */
78 +                free(ip->da);
79 +                ip->da = NULL;
80 +        }
81 +        ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
82 +        if (ip == NULL)
83 +                return(NULL);
84 +        ip->ns = nsamps;
85 +        return(ip);
86 + }
87 +
88 + /* Set minimum distance under which samples will start to merge */
89 + void
90 + interp2_spacing(INTERP2 *ip, double mind)
91 + {
92 +        if (mind <= 0)
93 +                return;
94 +        if ((.998*ip->dmin <= mind) & (mind <= 1.002*ip->dmin))
95 +                return;
96 +        if (ip->da != NULL) {   /* will need to recompute distribution */
97 +                free(ip->da);
98 +                ip->da = NULL;
99 +        }
100 +        ip->dmin = mind;
101 + }
102 +
103 + /* Modify smoothing parameter by the given factor */
104 + void
105 + interp2_smooth(INTERP2 *ip, double sf)
106 + {
107 +        if ((ip->smf *= sf) < NI2DSMF)
108 +                ip->smf = NI2DSMF;
109 + }
110 +
111   /* private call-back to sort position index */
112   static int
113   cmp_spos(const void *p1, const void *p2)
# Line 74 | Line 122 | cmp_spos(const void *p1, const void *p2)
122          return 0;
123   }
124  
125 < /* private routine to encode radius with range checks */
125 > /* private routine to order samples in a particular direction */
126 > static void
127 > sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang)
128 > {
129 >        const double    cosd = cos(ang);
130 >        const double    sind = sin(ang);
131 >        int             i;
132 >
133 >        for (i = ip->ns; i--; ) {
134 >                sord[i].si = i;
135 >                sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
136 >        }
137 >        qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
138 > }
139 >
140 > /* private routine to encode sample diameter with range checks */
141   static int
142 < encode_radius(const INTERP2 *ip, double r)
142 > encode_diameter(const INTERP2 *ip, double d)
143   {
144 <        const int       er = ENCODE_RAD(ip, r);
144 >        const int       ed = ENCODE_DIA(ip, d);
145  
146 <        if (er <= 0)
146 >        if (ed <= 0)
147                  return(0);
148 <        if (er >= 0xffff)
148 >        if (ed >= 0xffff)
149                  return(0xffff);
150 <        return(er);
150 >        return(ed);
151   }
152  
153 < /* Compute anisotropic Gaussian basis function interpolant */
154 < static int
155 < interp2_compute(INTERP2 *ip)
153 > /* (Re)compute anisotropic basis function interpolant (normally automatic) */
154 > int
155 > interp2_analyze(INTERP2 *ip)
156   {
157          SAMPORD *sortord;
158 <        int     *rightrndx, *leftrndx;
159 <        int     bd;
158 >        int     *rightrndx, *leftrndx, *endrndx;
159 >        int     i, bd;
160                                          /* sanity checks */
161 <        if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
161 >        if (ip == NULL)
162                  return(0);
163 <                                        /* need to allocate? */
164 <        if (ip->ra == NULL) {
165 <                ip->ra = (unsigned short (*)[NI2DIR])malloc(
103 <                                sizeof(unsigned short)*NI2DIR*ip->ns);
104 <                if (ip->ra == NULL)
105 <                        return(0);
163 >        if (ip->da != NULL) {           /* free previous data if any */
164 >                free(ip->da);
165 >                ip->da = NULL;
166          }
167 +        if ((ip->ns <= 1) | (ip->dmin <= 0))
168 +                return(0);
169 +                                        /* compute sample domain */
170 +        ip->smin[0] = ip->smin[1] = FHUGE;
171 +        ip->smax[0] = ip->smax[1] = -FHUGE;
172 +        for (i = ip->ns; i--; ) {
173 +                if (ip->spt[i][0] < ip->smin[0])
174 +                        ip->smin[0] = ip->spt[i][0];
175 +                if (ip->spt[i][0] > ip->smax[0])
176 +                        ip->smax[0] = ip->spt[i][0];
177 +                if (ip->spt[i][1] < ip->smin[1])
178 +                        ip->smin[1] = ip->spt[i][1];
179 +                if (ip->spt[i][1] > ip->smax[1])
180 +                        ip->smax[1] = ip->spt[i][1];
181 +        }
182 +        ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
183 +                        (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
184 +                                (1./NI2DIM/NI2DIM);
185 +        if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
186 +                return(0);
187 +                                        /* allocate analysis data */
188 +        ip->da = (struct interp2_samp *)calloc( ip->ns,
189 +                                        sizeof(struct interp2_samp) );
190 +        if (ip->da == NULL)
191 +                return(0);
192                                          /* get temporary arrays */
193          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
194          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
195          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
196 <        if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL))
196 >        endrndx = (int *)malloc(sizeof(int)*ip->ns);
197 >        if ((sortord == NULL) | (rightrndx == NULL) |
198 >                        (leftrndx == NULL) | (endrndx == NULL))
199                  return(0);
200                                          /* run through bidirections */
201          for (bd = 0; bd < NI2DIR/2; bd++) {
202              const double        ang = 2.*PI/NI2DIR*bd;
203 <            double              cosd, sind;
117 <            int                 i;
203 >            int                 *sptr;
204                                          /* create right reverse index */
205 <            if (bd) {                   /* re-use from prev. iteration? */
206 <                int     *sptr = rightrndx;
205 >            if (bd) {                   /* re-use from previous iteration? */
206 >                sptr = rightrndx;
207                  rightrndx = leftrndx;
208                  leftrndx = sptr;
209 <            } else {                    /* else compute it */
210 <                cosd = cos(ang + (PI/2. - PI/NI2DIR));
211 <                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++)
209 >            } else {                    /* else sort first half-plane */
210 >                sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
211 >                for (i = ip->ns; i--; )
212                      rightrndx[sortord[i].si] = i;
213 +                                        /* & store reverse order for later */
214 +                for (i = ip->ns; i--; )
215 +                    endrndx[sortord[i].si] = ip->ns-1 - i;
216              }
217                                          /* create new left reverse index */
218 <            cosd = cos(ang + (PI/2. + PI/NI2DIR));
219 <            sind = sin(ang + (PI/2. + PI/NI2DIR));
220 <            for (i = 0; i < ip->ns; i++) {
221 <                    sortord[i].si = i;
222 <                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
223 <            }
224 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
142 <            for (i = 0; i < ip->ns; i++)
218 >            if (bd == NI2DIR/2 - 1) {   /* use order from first iteration? */
219 >                sptr = leftrndx;
220 >                leftrndx = endrndx;
221 >                endrndx = sptr;
222 >            } else {                    /* else compute new half-plane */
223 >                sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
224 >                for (i = ip->ns; i--; )
225                      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];
226              }
227 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
227 >                                        /* sort grid values in this direction */
228 >            sort_samples(sortord, ip, ang);
229                                          /* find nearest neighbors each side */
230 <            for (i = 0; i < ip->ns; i++) {
231 <                const int       rpos = rightrndx[sortord[i].si];
155 <                const int       lpos = leftrndx[sortord[i].si];
230 >            for (i = ip->ns; i--; ) {
231 >                const int       ii = sortord[i].si;
232                  int             j;
233 <                                        /* preload with large radius */
234 <                ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
235 <                            .25*(sortord[ip->ns-1].dm - sortord[0].dm));
233 >                                        /* preload with large radii */
234 >                ip->da[ii].dia[bd] =
235 >                ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
236 >                                .5*(sortord[ip->ns-1].dm - sortord[0].dm));
237                  for (j = i; ++j < ip->ns; )     /* nearest above */
238 <                    if (rightrndx[sortord[j].si] > rpos &&
239 <                                    leftrndx[sortord[j].si] < lpos) {
240 <                        ip->ra[i][bd] = encode_radius(ip,
241 <                                        .5*(sortord[j].dm - sortord[i].dm));
238 >                    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
239 >                                    leftrndx[sortord[j].si] < leftrndx[ii]) {
240 >                        ip->da[ii].dia[bd] = encode_diameter(ip,
241 >                                                sortord[j].dm - sortord[i].dm);
242                          break;
243                      }
244                  for (j = i; j-- > 0; )          /* nearest below */
245 <                    if (rightrndx[sortord[j].si] < rpos &&
246 <                                    leftrndx[sortord[j].si] > lpos) {
247 <                        ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
248 <                                        .5*(sortord[i].dm - sortord[j].dm));
245 >                    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
246 >                                    leftrndx[sortord[j].si] > leftrndx[ii]) {
247 >                        ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
248 >                                                sortord[i].dm - sortord[j].dm);
249                          break;
250                      }
251              }
# Line 176 | Line 253 | interp2_compute(INTERP2 *ip)
253          free(sortord);                  /* clean up */
254          free(rightrndx);
255          free(leftrndx);
256 +        free(endrndx);
257          return(1);
258   }
259  
260 < /* private call returns log of raw weight for a particular sample */
261 < static double
262 < get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
260 > /* Compute unnormalized weight for a position relative to a sample */
261 > double
262 > interp2_wti(INTERP2 *ip, const int i, double x, double y)
263   {
264 <        double  dir, rd;
265 <        int     ri;
264 >        double          dir, rd, r2, d2;
265 >        int             ri;
266                                  /* get relative direction */
267          x -= ip->spt[i][0];
268          y -= ip->spt[i][1];
# Line 194 | Line 272 | get_ln_wt(const INTERP2 *ip, const int i, double x, do
272          rd = dir * (NI2DIR/2./PI);
273          ri = (int)rd;
274          rd -= (double)ri;
275 <        rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
276 <        rd = ip->smf * DECODE_RAD(ip, rd);
277 <                                /* return log of Gaussian weight */
278 <        return( (x*x + y*y) / (-2.*rd*rd) );
275 >        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
276 >        rd = ip->smf * DECODE_DIA(ip, rd);
277 >        r2 = 2.*rd*rd;
278 >        d2 = x*x + y*y;
279 >        if (d2 > 21.*r2)        /* result would be < 1e-9 */
280 >                return(.0);
281 >                                /* Gaussian times harmonic weighting */
282 >        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
283   }
284  
285 + /* private call to get grid flag index */
286 + static int
287 + interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
288 + {
289 +        int     ingrid = 1;
290 +
291 +        if (ip == NULL)         /* paranoia */
292 +                return(-1);
293 +                                /* need to compute interpolant? */
294 +        if (ip->da == NULL && !interp2_analyze(ip))
295 +                return(-1);
296 +                                /* get grid position */
297 +        fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
298 +        if (fgi[0] >= NI2DIM) {
299 +                fgi[0] = NI2DIM-1;
300 +                ingrid = 0;
301 +        } else if (fgi[0] < 0) {
302 +                fgi[0] = 0;
303 +                ingrid = 0;
304 +        }
305 +        fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
306 +        if (fgi[1] >= NI2DIM) {
307 +                fgi[1] = NI2DIM-1;
308 +                ingrid = 0;
309 +        } else if (fgi[1] < 0) {
310 +                fgi[1] = 0;
311 +                ingrid = 0;
312 +        }
313 +        return(ingrid);
314 + }
315 +
316 + /* private call to set black flag if not too close to the given sample */
317 + static void
318 + setblk(INTERP2 *ip, const int i, const int gi[2])
319 + {
320 +        double  dx = (gi[0]+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
321 +                                        ip->smin[0] - ip->spt[i][0];
322 +        double  dy = (gi[1]+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
323 +                                        ip->smin[1] - ip->spt[i][1];
324 +
325 +        if (dx*dx + dy*dy > 2.*ip->grid2)
326 +                ip->da[i].blkflg[gi[1]] |= 1<<gi[0];
327 + }
328 +
329 + #define chkblk(ip,i,gi) ((ip)->da[i].blkflg[(gi)[1]]>>(gi)[0] & 1)
330 +
331   /* Assign full set of normalized weights to interpolate the given position */
332   int
333   interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
334   {
335          double  wnorm;
336 +        int     fgi[2];
337 +        int     ingrid;
338          int     i;
339  
340 <        if ((wtv == NULL) | (ip == NULL))
340 >        if (wtv == NULL)
341                  return(0);
342 <                                        /* need to compute interpolant? */
343 <        if (ip->ra == NULL && !interp2_compute(ip))
342 >                                        /* get flag position */
343 >        if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
344                  return(0);
345  
346          wnorm = 0;                      /* compute raw weights */
347 <        for (i = ip->ns; i--; ) {
348 <                double  wt = get_ln_wt(ip, i, x, y);
349 <                if (wt < -21.) {
350 <                        wtv[i] = 0;     /* ignore weights < 1e-9 */
351 <                        continue;
222 <                }
223 <                wt = exp(wt);           /* Gaussian weight */
347 >        for (i = ip->ns; i--; )
348 >            if (chkblk(ip, i, fgi)) {
349 >                wtv[i] = 0;
350 >            } else {
351 >                double  wt = interp2_wti(ip, i, x, y);
352                  wtv[i] = wt;
353                  wnorm += wt;
354 <        }
354 >                if (wt <= 1e-9 && ingrid)
355 >                        setblk(ip, i, fgi);
356 >            }
357          if (wnorm <= 0)                 /* too far from all our samples! */
358                  return(0);
359          wnorm = 1./wnorm;               /* normalize weights */
# Line 238 | Line 368 | int
368   interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
369   {
370          int     nn = 0;
371 +        int     fgi[2];
372 +        int     ingrid;
373          double  wnorm;
374          int     i, j;
375  
376 <        if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
376 >        if ((n <= 0) | (wt == NULL) | (si == NULL))
377                  return(0);
378 <                                        /* need to compute interpolant? */
379 <        if (ip->ra == NULL && !interp2_compute(ip))
378 >                                        /* get flag position */
379 >        if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
380                  return(0);
381                                          /* identify top n weights */
382 <        for (i = ip->ns; i--; ) {
383 <                const double    lnwt = get_ln_wt(ip, i, x, y);
382 >        for (i = ip->ns; i--; )
383 >            if (!chkblk(ip, i, fgi)) {
384 >                const double    wti = interp2_wti(ip, i, x, y);
385 >                if (wti <= 1e-9) {
386 >                        if (ingrid)
387 >                                setblk(ip, i, fgi);
388 >                        continue;
389 >                }
390                  for (j = nn; j > 0; j--) {
391 <                        if (wt[j-1] >= lnwt)
391 >                        if (wt[j-1] >= wti)
392                                  break;
393                          if (j < n) {
394                                  wt[j] = wt[j-1];
# Line 258 | Line 396 | interp2_topsamp(float wt[], int si[], const int n, INT
396                          }
397                  }
398                  if (j < n) {            /* add/insert sample */
399 <                        wt[j] = lnwt;
399 >                        wt[j] = wti;
400                          si[j] = i;
401                          nn += (nn < n);
402                  }
403 <        }
404 <        wnorm = 0;                      /* exponentiate and normalize */
405 <        for (j = nn; j--; ) {
406 <                double  dwt = exp(wt[j]);
269 <                wt[j] = dwt;
270 <                wnorm += dwt;
271 <        }
403 >            }
404 >        wnorm = 0;                      /* normalize sample weights */
405 >        for (j = nn; j--; )
406 >                wnorm += wt[j];
407          if (wnorm <= 0)
408                  return(0);
409          wnorm = 1./wnorm;
# Line 283 | Line 418 | interp2_free(INTERP2 *ip)
418   {
419          if (ip == NULL)
420                  return;
421 <        if (ip->ra != NULL)
422 <                free(ip->ra);
421 >        if (ip->da != NULL)
422 >                free(ip->da);
423          free(ip);
424   }

Diff Legend

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