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.6 by greg, Tue Feb 12 02:56:05 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 39 | Line 40 | typedef struct {
40          float   dm;             /* distance measure in this direction */
41   } SAMPORD;
42  
43 < /* Allocate a new set of interpolation samples */
43 > /* Allocate a new set of interpolation samples (caller assigns spt[] array) */
44   INTERP2 *
45   interp2_alloc(int nsamps)
46   {
# 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   }
63  
64 + /* 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 +        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)
82 +                return(NULL);
83 +        ip->ns = nsamps;
84 +        return(ip);
85 + }
86 +
87 + /* 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 ((.998*ip->dmin <= mind) && (mind <= 1.002*ip->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   /* private call-back to sort position index */
111   static int
112   cmp_spos(const void *p1, const void *p2)
# Line 74 | Line 121 | cmp_spos(const void *p1, const void *p2)
121          return 0;
122   }
123  
124 < /* private routine to encode radius with range checks */
124 > /* 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 = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
135 >        }
136 >        qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
137 > }
138 >
139 > /* private routine to encode sample diameter with range checks */
140   static int
141 < encode_radius(const INTERP2 *ip, double r)
141 > encode_diameter(const INTERP2 *ip, double d)
142   {
143 <        const int       er = ENCODE_RAD(ip, r);
143 >        const int       ed = ENCODE_DIA(ip, d);
144  
145 <        if (er <= 0)
145 >        if (ed <= 0)
146                  return(0);
147 <        if (er >= 0xffff)
147 >        if (ed >= 0xffff)
148                  return(0xffff);
149 <        return(er);
149 >        return(ed);
150   }
151  
152 < /* Compute anisotropic Gaussian basis function interpolant */
153 < static int
154 < interp2_compute(INTERP2 *ip)
152 > /* (Re)compute anisotropic basis function interpolant (normally automatic) */
153 > int
154 > interp2_analyze(INTERP2 *ip)
155   {
156          SAMPORD *sortord;
157 <        int     *rightrndx, *leftrndx;
157 >        int     *rightrndx, *leftrndx, *endrndx;
158          int     bd;
159                                          /* sanity checks */
160 <        if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
160 >        if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
161                  return(0);
162                                          /* need to allocate? */
163 <        if (ip->ra == NULL) {
164 <                ip->ra = (unsigned short (*)[NI2DIR])malloc(
163 >        if (ip->da == NULL) {
164 >                ip->da = (unsigned short (*)[NI2DIR])malloc(
165                                  sizeof(unsigned short)*NI2DIR*ip->ns);
166 <                if (ip->ra == NULL)
166 >                if (ip->da == NULL)
167                          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 <        if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL))
173 >        endrndx = (int *)malloc(sizeof(int)*ip->ns);
174 >        if ((sortord == NULL) | (rightrndx == NULL) |
175 >                        (leftrndx == NULL) | (endrndx == NULL))
176                  return(0);
177                                          /* run through bidirections */
178          for (bd = 0; bd < NI2DIR/2; bd++) {
179              const double        ang = 2.*PI/NI2DIR*bd;
180 <            double              cosd, sind;
180 >            int                 *sptr;
181              int                 i;
182                                          /* create right reverse index */
183 <            if (bd) {                   /* re-use from prev. iteration? */
184 <                int     *sptr = rightrndx;
183 >            if (bd) {                   /* re-use from previous iteration? */
184 >                sptr = rightrndx;
185                  rightrndx = leftrndx;
186                  leftrndx = sptr;
187 <            } else {                    /* else compute it */
188 <                cosd = cos(ang + (PI/2. - PI/NI2DIR));
189 <                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++)
187 >            } else {                    /* else sort first half-plane */
188 >                sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
189 >                for (i = ip->ns; i--; )
190                      rightrndx[sortord[i].si] = i;
191 +                                        /* & store reverse order for later */
192 +                for (i = ip->ns; i--; )
193 +                    endrndx[sortord[i].si] = ip->ns-1 - i;
194              }
195                                          /* create new left reverse index */
196 <            cosd = cos(ang + (PI/2. + PI/NI2DIR));
197 <            sind = sin(ang + (PI/2. + PI/NI2DIR));
198 <            for (i = 0; i < ip->ns; i++) {
199 <                    sortord[i].si = i;
200 <                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
201 <            }
202 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
142 <            for (i = 0; i < ip->ns; i++)
196 >            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;
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];
204              }
205 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
205 >                                        /* sort grid values in this direction */
206 >            sort_samples(sortord, ip, ang);
207                                          /* find nearest neighbors each side */
208 <            for (i = 0; i < ip->ns; i++) {
209 <                const int       rpos = rightrndx[sortord[i].si];
155 <                const int       lpos = leftrndx[sortord[i].si];
208 >            for (i = ip->ns; i--; ) {
209 >                const int       ii = sortord[i].si;
210                  int             j;
211 <                                        /* preload with large radius */
212 <                ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
213 <                            .25*(sortord[ip->ns-1].dm - sortord[0].dm));
211 >                                        /* preload with large radii */
212 >                ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
213 >                            .5*(sortord[ip->ns-1].dm - sortord[0].dm));
214                  for (j = i; ++j < ip->ns; )     /* nearest above */
215 <                    if (rightrndx[sortord[j].si] > rpos &&
216 <                                    leftrndx[sortord[j].si] < lpos) {
217 <                        ip->ra[i][bd] = encode_radius(ip,
218 <                                        .5*(sortord[j].dm - sortord[i].dm));
215 >                    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
216 >                                    leftrndx[sortord[j].si] < leftrndx[ii]) {
217 >                        ip->da[ii][bd] = encode_diameter(ip,
218 >                                                sortord[j].dm - sortord[i].dm);
219                          break;
220                      }
221                  for (j = i; j-- > 0; )          /* nearest below */
222 <                    if (rightrndx[sortord[j].si] < rpos &&
223 <                                    leftrndx[sortord[j].si] > lpos) {
224 <                        ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
225 <                                        .5*(sortord[i].dm - sortord[j].dm));
222 >                    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
223 >                                    leftrndx[sortord[j].si] > leftrndx[ii]) {
224 >                        ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
225 >                                                sortord[i].dm - sortord[j].dm);
226                          break;
227                      }
228              }
# Line 176 | Line 230 | interp2_compute(INTERP2 *ip)
230          free(sortord);                  /* clean up */
231          free(rightrndx);
232          free(leftrndx);
233 +        free(endrndx);
234          return(1);
235   }
236  
237 < /* private call returns log of raw weight for a particular sample */
237 > /* private call returns raw weight for a particular sample */
238   static double
239 < get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
239 > get_wt(const INTERP2 *ip, const int i, double x, double y)
240   {
241 <        double  dir, rd;
241 >        double  dir, rd, r2, d2;
242          int     ri;
243                                  /* get relative direction */
244          x -= ip->spt[i][0];
# Line 194 | Line 249 | get_ln_wt(const INTERP2 *ip, const int i, double x, do
249          rd = dir * (NI2DIR/2./PI);
250          ri = (int)rd;
251          rd -= (double)ri;
252 <        rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
253 <        rd = ip->smf * DECODE_RAD(ip, rd);
254 <                                /* return log of Gaussian weight */
255 <        return( (x*x + y*y) / (-2.*rd*rd) );
252 >        rd = (1.-rd)*ip->da[i][ri] + rd*ip->da[i][(ri+1)%NI2DIR];
253 >        rd = ip->smf * DECODE_DIA(ip, rd);
254 >        r2 = 2.*rd*rd;
255 >        d2 = x*x + y*y;
256 >        if (d2 > 21.*r2)        /* result would be < 1e-9 */
257 >                return(.0);
258 >                                /* Gaussian times harmonic weighting */
259 >        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
260   }
261  
262   /* Assign full set of normalized weights to interpolate the given position */
# Line 210 | Line 269 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
269          if ((wtv == NULL) | (ip == NULL))
270                  return(0);
271                                          /* need to compute interpolant? */
272 <        if (ip->ra == NULL && !interp2_compute(ip))
272 >        if (ip->da == NULL && !interp2_analyze(ip))
273                  return(0);
274  
275          wnorm = 0;                      /* compute raw weights */
276          for (i = ip->ns; i--; ) {
277 <                double  wt = get_ln_wt(ip, i, x, y);
219 <                if (wt < -21.) {
220 <                        wtv[i] = 0;     /* ignore weights < 1e-9 */
221 <                        continue;
222 <                }
223 <                wt = exp(wt);           /* Gaussian weight */
277 >                double  wt = get_wt(ip, i, x, y);
278                  wtv[i] = wt;
279                  wnorm += wt;
280          }
# Line 244 | Line 298 | interp2_topsamp(float wt[], int si[], const int n, INT
298          if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
299                  return(0);
300                                          /* need to compute interpolant? */
301 <        if (ip->ra == NULL && !interp2_compute(ip))
301 >        if (ip->da == NULL && !interp2_analyze(ip))
302                  return(0);
303                                          /* identify top n weights */
304          for (i = ip->ns; i--; ) {
305 <                const double    lnwt = get_ln_wt(ip, i, x, y);
305 >                const double    wti = get_wt(ip, i, x, y);
306                  for (j = nn; j > 0; j--) {
307 <                        if (wt[j-1] >= lnwt)
307 >                        if (wt[j-1] >= wti)
308                                  break;
309                          if (j < n) {
310                                  wt[j] = wt[j-1];
# Line 258 | Line 312 | interp2_topsamp(float wt[], int si[], const int n, INT
312                          }
313                  }
314                  if (j < n) {            /* add/insert sample */
315 <                        wt[j] = lnwt;
315 >                        wt[j] = wti;
316                          si[j] = i;
317                          nn += (nn < n);
318                  }
319          }
320 <        wnorm = 0;                      /* exponentiate and normalize */
321 <        for (j = nn; j--; ) {
322 <                double  dwt = exp(wt[j]);
269 <                wt[j] = dwt;
270 <                wnorm += dwt;
271 <        }
320 >        wnorm = 0;                      /* normalize sample weights */
321 >        for (j = nn; j--; )
322 >                wnorm += wt[j];
323          if (wnorm <= 0)
324                  return(0);
325          wnorm = 1./wnorm;
# Line 283 | Line 334 | interp2_free(INTERP2 *ip)
334   {
335          if (ip == NULL)
336                  return;
337 <        if (ip->ra != NULL)
338 <                free(ip->ra);
337 >        if (ip->da != NULL)
338 >                free(ip->da);
339          free(ip);
340   }

Diff Legend

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