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.10 by greg, Thu Feb 14 19:57:10 2013 UTC vs.
Revision 2.17 by greg, Fri Mar 19 21:16:15 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 and an array
27 < * of flags indicating where weights are (nearly) zero.
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 41 | 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 57 | 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 66 | 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) {
# Line 79 | Line 95 | interp2_realloc(INTERP2 *ip, int nsamps)
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 >                        free(old_ip);
103 >                        return(NULL);
104 >                }
105 >        }
106          ip->ns = nsamps;
107          return(ip);
108   }
# Line 100 | Line 122 | interp2_spacing(INTERP2 *ip, double mind)
122          ip->dmin = mind;
123   }
124  
125 + /* Compute unnormalized weight for a position relative to a sample */
126 + double
127 + interp2_wti(INTERP2 *ip, const int i, double x, double y)
128 + {
129 +        double          dir, rd, r2, d2;
130 +        int             ri;
131 +                                /* get relative direction */
132 +        x -= ip->spt[i][0];
133 +        y -= ip->spt[i][1];
134 +        dir = atan2a(y, x);
135 +        dir += 2.*PI*(dir < 0);
136 +                                /* linear radius interpolation */
137 +        rd = dir * (NI2DIR/2./PI);
138 +        ri = (int)rd;
139 +        rd -= (double)ri;
140 +        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
141 +        rd = ip->smf * DECODE_DIA(ip, rd);
142 +        r2 = 2.*rd*rd;
143 +        d2 = x*x + y*y;
144 +        if (d2 > 21.*r2)        /* result would be < 1e-9 */
145 +                return(.0);
146 +                                /* Gaussian times harmonic weighting */
147 +        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
148 + }
149 +
150 + /* private call to get grid flag index */
151 + static int
152 + interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
153 + {
154 +        int     inbounds = 0;
155 +
156 +        if (ip == NULL)         /* paranoia */
157 +                return(-1);
158 +                                /* need to compute interpolant? */
159 +        if (ip->da == NULL && !interp2_analyze(ip))
160 +                return(-1);
161 +                                /* get x & y grid positions */
162 +        fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
163 +
164 +        if (fgi[0] >= NI2DIM)
165 +                fgi[0] = NI2DIM-1;
166 +        else if (fgi[0] < 0)
167 +                fgi[0] = 0;
168 +        else
169 +                ++inbounds;
170 +
171 +        fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
172 +
173 +        if (fgi[1] >= NI2DIM)
174 +                fgi[1] = NI2DIM-1;
175 +        else if (fgi[1] < 0)
176 +                fgi[1] = 0;
177 +        else
178 +                ++inbounds;
179 +
180 +        return(inbounds == 2);
181 + }
182 +
183 + #define setflg(fl,xi,yi)        ((fl)[yi] |= 1<<(xi))
184 +
185 + #define chkflg(fl,xi,yi)        ((fl)[yi]>>(xi) & 1)
186 +
187 + /* private flood function to determine sample influence */
188 + static void
189 + influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
190 +                        int xfi, int yfi)
191 + {
192 +        double  gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
193 +                                        ip->smin[0];
194 +        double  gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
195 +                                        ip->smin[1];
196 +        double  dx = gx - ip->spt[i][0];
197 +        double  dy = gy - ip->spt[i][1];
198 +
199 +        setflg(visited, xfi, yfi);
200 +
201 +        if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
202 +                return;
203 +
204 +        setflg(ip->da[i].infl, xfi, yfi);
205 +
206 +        if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
207 +                influence_flood(ip, i, visited, xfi-1, yfi);
208 +
209 +        if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
210 +                influence_flood(ip, i, visited, xfi, yfi-1);
211 +
212 +        if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
213 +                influence_flood(ip, i, visited, xfi+1, yfi);
214 +
215 +        if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
216 +                influence_flood(ip, i, visited, xfi, yfi+1);
217 + }
218 +
219 + /* private call to compute sample influence maps */
220 + static void
221 + map_influence(INTERP2 *ip)
222 + {
223 +        unsigned short  visited[NI2DIM];
224 +        int             fgi[2];
225 +        int             i, j;
226 +
227 +        for (i = ip->ns; i--; ) {
228 +                for (j = NI2DIM; j--; ) {
229 +                        ip->da[i].infl[j] = 0;
230 +                        visited[j] = 0;
231 +                }
232 +                interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
233 +
234 +                influence_flood(ip, i, visited, fgi[0], fgi[1]);
235 +        }
236 + }
237 +
238   /* Modify smoothing parameter by the given factor */
239   void
240   interp2_smooth(INTERP2 *ip, double sf)
241   {
242 +        float   old_smf = ip->smf;
243 +
244          if ((ip->smf *= sf) < NI2DSMF)
245                  ip->smf = NI2DSMF;
246 +                                        /* need to recompute influence maps? */
247 +        if (ip->da != NULL && (old_smf*.85 > ip->smf) |
248 +                                (ip->smf > old_smf*1.15))
249 +                map_influence(ip);
250   }
251  
252   /* private call-back to sort position index */
# Line 134 | Line 275 | sort_samples(SAMPORD *sord, const INTERP2 *ip, double
275                  sord[i].si = i;
276                  sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
277          }
278 <        qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
278 >        qsort(sord, ip->ns, sizeof(SAMPORD), cmp_spos);
279   }
280  
140 – /* private routine to encode sample diameter with range checks */
141 – static int
142 – encode_diameter(const INTERP2 *ip, double d)
143 – {
144 –        const int       ed = ENCODE_DIA(ip, d);
145 –
146 –        if (ed <= 0)
147 –                return(0);
148 –        if (ed >= 0xffff)
149 –                return(0xffff);
150 –        return(ed);
151 – }
152 –
281   /* (Re)compute anisotropic basis function interpolant (normally automatic) */
282   int
283   interp2_analyze(INTERP2 *ip)
# Line 168 | Line 296 | interp2_analyze(INTERP2 *ip)
296                  return(0);
297                                          /* compute sample domain */
298          ip->smin[0] = ip->smin[1] = FHUGE;
299 <        ip->smul[0] = ip->smul[1] = -FHUGE;
299 >        ip->smax[0] = ip->smax[1] = -FHUGE;
300          for (i = ip->ns; i--; ) {
301 <                if (ip->spt[i][0] < ip->smin[0])
302 <                        ip->smin[0] = ip->spt[i][0];
303 <                if (ip->spt[i][0] > ip->smul[0])
304 <                        ip->smul[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->smul[1])
180 <                        ip->smul[1] = ip->spt[i][1];
301 >                if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
302 >                if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
303 >                if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
304 >                if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
305          }
306 <        ip->smul[0] -= ip->smin[0];
307 <        ip->smul[1] -= ip->smin[1];
308 <        ip->grid2 = (ip->smul[0]*ip->smul[0] + ip->smul[1]*ip->smul[1]) *
185 <                        (4./NI2DIM/NI2DIM);
306 >        ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
307 >                        (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
308 >                                (1./NI2DIM/NI2DIM);
309          if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
310                  return(0);
188 –        if (ip->smul[0] > FTINY)
189 –                ip->smul[0] = NI2DIM / ip->smul[0];
190 –        if (ip->smul[1] > FTINY)
191 –                ip->smul[1] = NI2DIM / ip->smul[1];
311                                          /* allocate analysis data */
312 <        ip->da = (struct interp2_samp *)calloc( ip->ns,
313 <                                        sizeof(struct interp2_samp) );
312 >        ip->da = (struct interp2_samp *)malloc(
313 >                                sizeof(struct interp2_samp)*ip->ns );
314          if (ip->da == NULL)
315                  return(0);
316 <                                        /* get temporary arrays */
316 >                                        /* allocate temporary arrays */
317          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
318          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
319          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
# Line 255 | Line 374 | interp2_analyze(INTERP2 *ip)
374                      }
375              }
376          }
377 <        free(sortord);                  /* clean up */
377 >        free(sortord);                  /* release temp arrays */
378          free(rightrndx);
379          free(leftrndx);
380          free(endrndx);
381 <        return(1);
381 >                                        /* map sample influence areas */
382 >        map_influence(ip);
383 >        return(1);                      /* all done */
384   }
385  
265 – /* Compute unnormalized weight for a position relative to a sample */
266 – double
267 – interp2_wti(INTERP2 *ip, const int i, double x, double y)
268 – {
269 –        int     xfi, yfi;
270 –        double  dir, rd, r2, d2;
271 –        int     ri;
272 –                                /* need to compute interpolant? */
273 –        if (ip->da == NULL && !interp2_analyze(ip))
274 –                return(0);
275 –                                /* get grid position */
276 –        xfi = (x - ip->smin[0]) * ip->smul[0];
277 –        if (xfi >= NI2DIM)
278 –                xfi = NI2DIM-1;
279 –        else
280 –                xfi *= (xfi >= 0);
281 –        yfi = (y - ip->smin[1]) * ip->smul[1];
282 –        if (yfi >= NI2DIM)
283 –                yfi = NI2DIM-1;
284 –        else
285 –                yfi *= (yfi >= 0);
286 –        x -= ip->spt[i][0];     /* check distance */
287 –        y -= ip->spt[i][1];
288 –        d2 = x*x + y*y;
289 –                                /* zero weight this zone? */
290 –        if (d2 > ip->grid2 && ip->da[i].blkflg[yfi] & 1<<xfi)
291 –                return(.0);
292 –
293 –        dir = atan2a(y, x);     /* get relative direction */
294 –        dir += 2.*PI*(dir < 0);
295 –                                /* linear radius interpolation */
296 –        rd = dir * (NI2DIR/2./PI);
297 –        ri = (int)rd;
298 –        rd -= (double)ri;
299 –        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
300 –        rd = ip->smf * DECODE_DIA(ip, rd);
301 –        r2 = 2.*rd*rd;
302 –        if (d2 > 21.*r2) {      /* result would be < 1e-9 */
303 –                ip->da[i].blkflg[yfi] |= 1<<xfi;
304 –                return(.0);
305 –        }
306 –                                /* Gaussian times harmonic weighting */
307 –        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
308 – }
309 –
386   /* Assign full set of normalized weights to interpolate the given position */
387   int
388   interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
389   {
390          double  wnorm;
391 +        int     fgi[2];
392          int     i;
393  
394 <        if ((wtv == NULL) | (ip == NULL))
394 >        if (wtv == NULL)
395                  return(0);
396 +                                        /* get flag position */
397 +        if (interp2_flagpos(fgi, ip, x, y) < 0)
398 +                return(0);
399  
400          wnorm = 0;                      /* compute raw weights */
401 <        for (i = ip->ns; i--; ) {
401 >        for (i = ip->ns; i--; )
402 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
403                  double  wt = interp2_wti(ip, i, x, y);
404                  wtv[i] = wt;
405                  wnorm += wt;
406 <        }
406 >            } else
407 >                wtv[i] = 0;
408          if (wnorm <= 0)                 /* too far from all our samples! */
409                  return(0);
410          wnorm = 1./wnorm;               /* normalize weights */
# Line 337 | Line 419 | int
419   interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
420   {
421          int     nn = 0;
422 +        int     fgi[2];
423          double  wnorm;
424          int     i, j;
425  
426 <        if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
426 >        if ((n <= 0) | (wt == NULL) | (si == NULL))
427                  return(0);
428 +                                        /* get flag position */
429 +        if (interp2_flagpos(fgi, ip, x, y) < 0)
430 +                return(0);
431                                          /* identify top n weights */
432 <        for (i = ip->ns; i--; ) {
432 >        for (i = ip->ns; i--; )
433 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
434                  const double    wti = interp2_wti(ip, i, x, y);
348 –                if (wti <= 1e-9)
349 –                        continue;
435                  for (j = nn; j > 0; j--) {
436                          if (wt[j-1] >= wti)
437                                  break;
# Line 360 | Line 445 | interp2_topsamp(float wt[], int si[], const int n, INT
445                          si[j] = i;
446                          nn += (nn < n);
447                  }
448 <        }
448 >            }
449          wnorm = 0;                      /* normalize sample weights */
450          for (j = nn; j--; )
451                  wnorm += wt[j];

Diff Legend

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