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.7 by greg, Tue Feb 12 17:47:58 2013 UTC vs.
Revision 2.16 by greg, Thu Mar 11 01:58:59 2021 UTC

# Line 23 | Line 23 | static const char RCSid[] = "$Id$";
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.
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>
# Line 40 | 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 56 | Line 70 | interp2_alloc(int nsamps)
70          nip->ns = nsamps;
71          nip->dmin = 1;          /* default minimum diameter */
72          nip->smf = NI2DSMF;     /* default smoothing factor */
73 +        nip->c_data = NULL;
74          nip->da = NULL;
75                                  /* caller must assign spt[] array */
76          return(nip);
# Line 65 | Line 80 | interp2_alloc(int nsamps)
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);
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 <                return(NULL);
98 >        if (ip == NULL) {
99 >                if (nsamps <= ip->ns)
100 >                        ip = old_ip;
101 >                else
102 >                        return(NULL);
103 >        }
104          ip->ns = nsamps;
105          return(ip);
106   }
# Line 99 | Line 120 | interp2_spacing(INTERP2 *ip, double mind)
120          ip->dmin = mind;
121   }
122  
123 + /* Compute unnormalized weight for a position relative to a sample */
124 + double
125 + interp2_wti(INTERP2 *ip, const int i, double x, double y)
126 + {
127 +        double          dir, rd, r2, d2;
128 +        int             ri;
129 +                                /* get relative direction */
130 +        x -= ip->spt[i][0];
131 +        y -= ip->spt[i][1];
132 +        dir = atan2a(y, x);
133 +        dir += 2.*PI*(dir < 0);
134 +                                /* linear radius interpolation */
135 +        rd = dir * (NI2DIR/2./PI);
136 +        ri = (int)rd;
137 +        rd -= (double)ri;
138 +        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
139 +        rd = ip->smf * DECODE_DIA(ip, rd);
140 +        r2 = 2.*rd*rd;
141 +        d2 = x*x + y*y;
142 +        if (d2 > 21.*r2)        /* result would be < 1e-9 */
143 +                return(.0);
144 +                                /* Gaussian times harmonic weighting */
145 +        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
146 + }
147 +
148 + /* private call to get grid flag index */
149 + static int
150 + interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
151 + {
152 +        int     inbounds = 0;
153 +
154 +        if (ip == NULL)         /* paranoia */
155 +                return(-1);
156 +                                /* need to compute interpolant? */
157 +        if (ip->da == NULL && !interp2_analyze(ip))
158 +                return(-1);
159 +                                /* get x & y grid positions */
160 +        fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
161 +
162 +        if (fgi[0] >= NI2DIM)
163 +                fgi[0] = NI2DIM-1;
164 +        else if (fgi[0] < 0)
165 +                fgi[0] = 0;
166 +        else
167 +                ++inbounds;
168 +
169 +        fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
170 +
171 +        if (fgi[1] >= NI2DIM)
172 +                fgi[1] = NI2DIM-1;
173 +        else if (fgi[1] < 0)
174 +                fgi[1] = 0;
175 +        else
176 +                ++inbounds;
177 +
178 +        return(inbounds == 2);
179 + }
180 +
181 + #define setflg(fl,xi,yi)        ((fl)[yi] |= 1<<(xi))
182 +
183 + #define chkflg(fl,xi,yi)        ((fl)[yi]>>(xi) & 1)
184 +
185 + /* private flood function to determine sample influence */
186 + static void
187 + influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
188 +                        int xfi, int yfi)
189 + {
190 +        double  gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
191 +                                        ip->smin[0];
192 +        double  gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
193 +                                        ip->smin[1];
194 +        double  dx = gx - ip->spt[i][0];
195 +        double  dy = gy - ip->spt[i][1];
196 +
197 +        setflg(visited, xfi, yfi);
198 +
199 +        if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
200 +                return;
201 +
202 +        setflg(ip->da[i].infl, xfi, yfi);
203 +
204 +        if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
205 +                influence_flood(ip, i, visited, xfi-1, yfi);
206 +
207 +        if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
208 +                influence_flood(ip, i, visited, xfi, yfi-1);
209 +
210 +        if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
211 +                influence_flood(ip, i, visited, xfi+1, yfi);
212 +
213 +        if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
214 +                influence_flood(ip, i, visited, xfi, yfi+1);
215 + }
216 +
217 + /* private call to compute sample influence maps */
218 + static void
219 + map_influence(INTERP2 *ip)
220 + {
221 +        unsigned short  visited[NI2DIM];
222 +        int             fgi[2];
223 +        int             i, j;
224 +
225 +        for (i = ip->ns; i--; ) {
226 +                for (j = NI2DIM; j--; ) {
227 +                        ip->da[i].infl[j] = 0;
228 +                        visited[j] = 0;
229 +                }
230 +                interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
231 +
232 +                influence_flood(ip, i, visited, fgi[0], fgi[1]);
233 +        }
234 + }
235 +
236   /* Modify smoothing parameter by the given factor */
237   void
238   interp2_smooth(INTERP2 *ip, double sf)
239   {
240 +        float   old_smf = ip->smf;
241 +
242          if ((ip->smf *= sf) < NI2DSMF)
243                  ip->smf = NI2DSMF;
244 +                                        /* need to recompute influence maps? */
245 +        if (ip->da != NULL && (old_smf*.85 > ip->smf) |
246 +                                (ip->smf > old_smf*1.15))
247 +                map_influence(ip);
248   }
249  
250   /* private call-back to sort position index */
# Line 133 | Line 273 | sort_samples(SAMPORD *sord, const INTERP2 *ip, double
273                  sord[i].si = i;
274                  sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
275          }
276 <        qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
276 >        qsort(sord, ip->ns, sizeof(SAMPORD), cmp_spos);
277   }
278  
139 – /* private routine to encode sample diameter with range checks */
140 – static int
141 – encode_diameter(const INTERP2 *ip, double d)
142 – {
143 –        const int       ed = ENCODE_DIA(ip, d);
144 –
145 –        if (ed <= 0)
146 –                return(0);
147 –        if (ed >= 0xffff)
148 –                return(0xffff);
149 –        return(ed);
150 – }
151 –
279   /* (Re)compute anisotropic basis function interpolant (normally automatic) */
280   int
281   interp2_analyze(INTERP2 *ip)
282   {
283          SAMPORD *sortord;
284          int     *rightrndx, *leftrndx, *endrndx;
285 <        int     bd;
285 >        int     i, bd;
286                                          /* sanity checks */
287 <        if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
287 >        if (ip == NULL)
288                  return(0);
289 <                                        /* need to allocate? */
290 <        if (ip->da == NULL) {
291 <                ip->da = (unsigned short (*)[NI2DIR])malloc(
165 <                                sizeof(unsigned short)*NI2DIR*ip->ns);
166 <                if (ip->da == NULL)
167 <                        return(0);
289 >        if (ip->da != NULL) {           /* free previous data if any */
290 >                free(ip->da);
291 >                ip->da = NULL;
292          }
293 <                                        /* get temporary arrays */
293 >        if ((ip->ns <= 1) | (ip->dmin <= 0))
294 >                return(0);
295 >                                        /* compute sample domain */
296 >        ip->smin[0] = ip->smin[1] = FHUGE;
297 >        ip->smax[0] = ip->smax[1] = -FHUGE;
298 >        for (i = ip->ns; i--; ) {
299 >                if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
300 >                if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
301 >                if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
302 >                if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
303 >        }
304 >        ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
305 >                        (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
306 >                                (1./NI2DIM/NI2DIM);
307 >        if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
308 >                return(0);
309 >                                        /* allocate analysis data */
310 >        ip->da = (struct interp2_samp *)malloc(
311 >                                sizeof(struct interp2_samp)*ip->ns );
312 >        if (ip->da == NULL)
313 >                return(0);
314 >                                        /* allocate temporary arrays */
315          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
316          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
317          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
# Line 178 | Line 323 | interp2_analyze(INTERP2 *ip)
323          for (bd = 0; bd < NI2DIR/2; bd++) {
324              const double        ang = 2.*PI/NI2DIR*bd;
325              int                 *sptr;
181 –            int                 i;
326                                          /* create right reverse index */
327              if (bd) {                   /* re-use from previous iteration? */
328                  sptr = rightrndx;
# Line 209 | Line 353 | interp2_analyze(INTERP2 *ip)
353                  const int       ii = sortord[i].si;
354                  int             j;
355                                          /* preload with large radii */
356 <                ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
357 <                            .5*(sortord[ip->ns-1].dm - sortord[0].dm));
356 >                ip->da[ii].dia[bd] =
357 >                ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
358 >                                .5*(sortord[ip->ns-1].dm - sortord[0].dm));
359                  for (j = i; ++j < ip->ns; )     /* nearest above */
360                      if (rightrndx[sortord[j].si] > rightrndx[ii] &&
361                                      leftrndx[sortord[j].si] < leftrndx[ii]) {
362 <                        ip->da[ii][bd] = encode_diameter(ip,
362 >                        ip->da[ii].dia[bd] = encode_diameter(ip,
363                                                  sortord[j].dm - sortord[i].dm);
364                          break;
365                      }
366                  for (j = i; j-- > 0; )          /* nearest below */
367                      if (rightrndx[sortord[j].si] < rightrndx[ii] &&
368                                      leftrndx[sortord[j].si] > leftrndx[ii]) {
369 <                        ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
369 >                        ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
370                                                  sortord[i].dm - sortord[j].dm);
371                          break;
372                      }
373              }
374          }
375 <        free(sortord);                  /* clean up */
375 >        free(sortord);                  /* release temp arrays */
376          free(rightrndx);
377          free(leftrndx);
378          free(endrndx);
379 <        return(1);
379 >                                        /* map sample influence areas */
380 >        map_influence(ip);
381 >        return(1);                      /* all done */
382   }
383  
237 – /* private call returns raw weight for a particular sample */
238 – static double
239 – get_wt(const INTERP2 *ip, const int i, double x, double y)
240 – {
241 –        double  dir, rd, r2, d2;
242 –        int     ri;
243 –                                /* get relative direction */
244 –        x -= ip->spt[i][0];
245 –        y -= ip->spt[i][1];
246 –        dir = atan2a(y, x);
247 –        dir += 2.*PI*(dir < 0);
248 –                                /* linear radius interpolation */
249 –        rd = dir * (NI2DIR/2./PI);
250 –        ri = (int)rd;
251 –        rd -= (double)ri;
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 –
384   /* Assign full set of normalized weights to interpolate the given position */
385   int
386   interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
387   {
388          double  wnorm;
389 +        int     fgi[2];
390          int     i;
391  
392 <        if ((wtv == NULL) | (ip == NULL))
392 >        if (wtv == NULL)
393                  return(0);
394 <                                        /* need to compute interpolant? */
395 <        if (ip->da == NULL && !interp2_analyze(ip))
394 >                                        /* get flag position */
395 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
396                  return(0);
397  
398          wnorm = 0;                      /* compute raw weights */
399 <        for (i = ip->ns; i--; ) {
400 <                double  wt = get_wt(ip, i, x, y);
399 >        for (i = ip->ns; i--; )
400 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
401 >                double  wt = interp2_wti(ip, i, x, y);
402                  wtv[i] = wt;
403                  wnorm += wt;
404 <        }
404 >            } else
405 >                wtv[i] = 0;
406          if (wnorm <= 0)                 /* too far from all our samples! */
407                  return(0);
408          wnorm = 1./wnorm;               /* normalize weights */
# Line 292 | Line 417 | int
417   interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
418   {
419          int     nn = 0;
420 +        int     fgi[2];
421          double  wnorm;
422          int     i, j;
423  
424 <        if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
424 >        if ((n <= 0) | (wt == NULL) | (si == NULL))
425                  return(0);
426 <                                        /* need to compute interpolant? */
427 <        if (ip->da == NULL && !interp2_analyze(ip))
426 >                                        /* get flag position */
427 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
428                  return(0);
429                                          /* identify top n weights */
430 <        for (i = ip->ns; i--; ) {
431 <                const double    wti = get_wt(ip, i, x, y);
430 >        for (i = ip->ns; i--; )
431 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
432 >                const double    wti = interp2_wti(ip, i, x, y);
433                  for (j = nn; j > 0; j--) {
434                          if (wt[j-1] >= wti)
435                                  break;
# Line 316 | Line 443 | interp2_topsamp(float wt[], int si[], const int n, INT
443                          si[j] = i;
444                          nn += (nn < n);
445                  }
446 <        }
446 >            }
447          wnorm = 0;                      /* normalize sample weights */
448          for (j = nn; j--; )
449                  wnorm += wt[j];

Diff Legend

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