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.4 by greg, Mon Feb 11 22:56:22 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 > * direction.  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.
27 > ****************************************************************/
28  
29   #include <stdio.h>
30   #include <stdlib.h>
31   #include "rtmath.h"
32   #include "interp2d.h"
33  
34 < #define DECODE_RAD(ip,er)       ((ip)->rmin*(1. + .5*(er)))
35 < #define ENCODE_RAD(ip,r)        ((int)(2.*(r)/(ip)->rmin) - 2)
34 > #define DECODE_DIA(ip,ed)       ((ip)->dmin*(1. + .5*(ed)))
35 > #define ENCODE_DIA(ip,d)        ((int)(2.*(d)/(ip)->dmin) - 2)
36  
37   /* Sample order (private) */
38   typedef struct {
# Line 53 | Line 54 | interp2_alloc(int nsamps)
54                  return(NULL);
55  
56          nip->ns = nsamps;
57 <        nip->rmin = .5;         /* default radius minimum */
57 >        nip->dmin = 1;          /* default minimum diameter */
58          nip->smf = NI2DSMF;     /* default smoothing factor */
59 <        nip->ra = NULL;
59 >        nip->da = NULL;
60                                  /* caller must assign spt[] array */
61          return(nip);
62   }
# Line 72 | Line 73 | interp2_realloc(INTERP2 *ip, int nsamps)
73          }
74          if (nsamps == ip->ns);
75                  return(ip);
76 <        if (ip->ra != NULL) {   /* will need to recompute distribution */
77 <                free(ip->ra);
78 <                ip->ra = NULL;
76 >        if (ip->da != NULL) {   /* will need to recompute distribution */
77 >                free(ip->da);
78 >                ip->da = NULL;
79          }
80          ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
81          if (ip == NULL)
# Line 112 | Line 113 | sort_samples(SAMPORD *sord, const INTERP2 *ip, double
113          qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
114   }
115  
116 < /* private routine to encode radius with range checks */
116 > /* private routine to encode sample diameter with range checks */
117   static int
118 < encode_radius(const INTERP2 *ip, double r)
118 > encode_diameter(const INTERP2 *ip, double d)
119   {
120 <        const int       er = ENCODE_RAD(ip, r);
120 >        const int       ed = ENCODE_DIA(ip, d);
121  
122 <        if (er <= 0)
122 >        if (ed <= 0)
123                  return(0);
124 <        if (er >= 0xffff)
124 >        if (ed >= 0xffff)
125                  return(0xffff);
126 <        return(er);
126 >        return(ed);
127   }
128  
129   /* (Re)compute anisotropic basis function interpolant (normally automatic) */
# Line 133 | Line 134 | interp2_analyze(INTERP2 *ip)
134          int     *rightrndx, *leftrndx, *endrndx;
135          int     bd;
136                                          /* sanity checks */
137 <        if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
137 >        if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
138                  return(0);
139                                          /* need to allocate? */
140 <        if (ip->ra == NULL) {
141 <                ip->ra = (unsigned short (*)[NI2DIR])malloc(
140 >        if (ip->da == NULL) {
141 >                ip->da = (unsigned short (*)[NI2DIR])malloc(
142                                  sizeof(unsigned short)*NI2DIR*ip->ns);
143 <                if (ip->ra == NULL)
143 >                if (ip->da == NULL)
144                          return(0);
145          }
146                                          /* get temporary arrays */
# Line 182 | Line 183 | interp2_analyze(INTERP2 *ip)
183              sort_samples(sortord, ip, ang);
184                                          /* find nearest neighbors each side */
185              for (i = ip->ns; i--; ) {
186 <                const int       rpos = rightrndx[sortord[i].si];
186 <                const int       lpos = leftrndx[sortord[i].si];
186 >                const int       ii = sortord[i].si;
187                  int             j;
188 <                                        /* preload with large radius */
189 <                ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
190 <                            .25*(sortord[ip->ns-1].dm - sortord[0].dm));
188 >                                        /* preload with large radii */
189 >                ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
190 >                            .5*(sortord[ip->ns-1].dm - sortord[0].dm));
191                  for (j = i; ++j < ip->ns; )     /* nearest above */
192 <                    if (rightrndx[sortord[j].si] > rpos &&
193 <                                    leftrndx[sortord[j].si] < lpos) {
194 <                        ip->ra[i][bd] = encode_radius(ip,
195 <                                        .5*(sortord[j].dm - sortord[i].dm));
192 >                    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
193 >                                    leftrndx[sortord[j].si] < leftrndx[ii]) {
194 >                        ip->da[ii][bd] = encode_diameter(ip,
195 >                                                sortord[j].dm - sortord[i].dm);
196                          break;
197                      }
198                  for (j = i; j-- > 0; )          /* nearest below */
199 <                    if (rightrndx[sortord[j].si] < rpos &&
200 <                                    leftrndx[sortord[j].si] > lpos) {
201 <                        ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
202 <                                        .5*(sortord[i].dm - sortord[j].dm));
199 >                    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
200 >                                    leftrndx[sortord[j].si] > leftrndx[ii]) {
201 >                        ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
202 >                                                sortord[i].dm - sortord[j].dm);
203                          break;
204                      }
205              }
# Line 211 | Line 211 | interp2_analyze(INTERP2 *ip)
211          return(1);
212   }
213  
214 < /* private call returns log of raw weight for a particular sample */
214 > /* private call returns raw weight for a particular sample */
215   static double
216 < get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
216 > get_wt(const INTERP2 *ip, const int i, double x, double y)
217   {
218 <        double  dir, rd;
218 >        double  dir, rd, d2;
219          int     ri;
220                                  /* get relative direction */
221          x -= ip->spt[i][0];
# Line 226 | Line 226 | get_ln_wt(const INTERP2 *ip, const int i, double x, do
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) );
229 >        rd = (1.-rd)*ip->da[i][ri] + rd*ip->da[i][(ri+1)%NI2DIR];
230 >        rd = ip->smf * DECODE_DIA(ip, rd);
231 >        d2 = x*x + y*y;
232 >                                /* Gaussian times harmonic weighting */
233 >        return( exp(d2/(-2.*rd*rd)) * ip->dmin/(ip->dmin + sqrt(d2)) );
234   }
235  
236   /* Assign full set of normalized weights to interpolate the given position */
# Line 242 | Line 243 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
243          if ((wtv == NULL) | (ip == NULL))
244                  return(0);
245                                          /* need to compute interpolant? */
246 <        if (ip->ra == NULL && !interp2_analyze(ip))
246 >        if (ip->da == NULL && !interp2_analyze(ip))
247                  return(0);
248  
249          wnorm = 0;                      /* compute raw weights */
250          for (i = ip->ns; i--; ) {
251 <                double  wt = get_ln_wt(ip, i, x, y);
251 <                if (wt < -21.) {
252 <                        wtv[i] = 0;     /* ignore weights < 1e-9 */
253 <                        continue;
254 <                }
255 <                wt = exp(wt);           /* Gaussian weight */
251 >                double  wt = get_wt(ip, i, x, y);
252                  wtv[i] = wt;
253                  wnorm += wt;
254          }
# Line 276 | Line 272 | interp2_topsamp(float wt[], int si[], const int n, INT
272          if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
273                  return(0);
274                                          /* need to compute interpolant? */
275 <        if (ip->ra == NULL && !interp2_analyze(ip))
275 >        if (ip->da == NULL && !interp2_analyze(ip))
276                  return(0);
277                                          /* identify top n weights */
278          for (i = ip->ns; i--; ) {
279 <                const double    lnwt = get_ln_wt(ip, i, x, y);
279 >                const double    wti = get_wt(ip, i, x, y);
280                  for (j = nn; j > 0; j--) {
281 <                        if (wt[j-1] >= lnwt)
281 >                        if (wt[j-1] >= wti)
282                                  break;
283                          if (j < n) {
284                                  wt[j] = wt[j-1];
# Line 290 | Line 286 | interp2_topsamp(float wt[], int si[], const int n, INT
286                          }
287                  }
288                  if (j < n) {            /* add/insert sample */
289 <                        wt[j] = lnwt;
289 >                        wt[j] = wti;
290                          si[j] = i;
291                          nn += (nn < n);
292                  }
293          }
294 <        wnorm = 0;                      /* exponentiate and normalize */
295 <        for (j = nn; j--; ) {
296 <                double  dwt = exp(wt[j]);
301 <                wt[j] = dwt;
302 <                wnorm += dwt;
303 <        }
294 >        wnorm = 0;                      /* normalize sample weights */
295 >        for (j = nn; j--; )
296 >                wnorm += wt[j];
297          if (wnorm <= 0)
298                  return(0);
299          wnorm = 1./wnorm;
# Line 315 | Line 308 | interp2_free(INTERP2 *ip)
308   {
309          if (ip == NULL)
310                  return;
311 <        if (ip->ra != NULL)
312 <                free(ip->ra);
311 >        if (ip->da != NULL)
312 >                free(ip->da);
313          free(ip);
314   }

Diff Legend

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