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.11 by greg, Fri Feb 15 01:26:47 2013 UTC vs.
Revision 2.15 by greg, Sat Feb 13 16:49:18 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 100 | Line 114 | interp2_spacing(INTERP2 *ip, double mind)
114          ip->dmin = mind;
115   }
116  
117 + /* Compute unnormalized weight for a position relative to a sample */
118 + double
119 + interp2_wti(INTERP2 *ip, const int i, double x, double y)
120 + {
121 +        double          dir, rd, r2, d2;
122 +        int             ri;
123 +                                /* get relative direction */
124 +        x -= ip->spt[i][0];
125 +        y -= ip->spt[i][1];
126 +        dir = atan2a(y, x);
127 +        dir += 2.*PI*(dir < 0);
128 +                                /* linear radius interpolation */
129 +        rd = dir * (NI2DIR/2./PI);
130 +        ri = (int)rd;
131 +        rd -= (double)ri;
132 +        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
133 +        rd = ip->smf * DECODE_DIA(ip, rd);
134 +        r2 = 2.*rd*rd;
135 +        d2 = x*x + y*y;
136 +        if (d2 > 21.*r2)        /* result would be < 1e-9 */
137 +                return(.0);
138 +                                /* Gaussian times harmonic weighting */
139 +        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
140 + }
141 +
142 + /* private call to get grid flag index */
143 + static int
144 + interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
145 + {
146 +        int     inbounds = 0;
147 +
148 +        if (ip == NULL)         /* paranoia */
149 +                return(-1);
150 +                                /* need to compute interpolant? */
151 +        if (ip->da == NULL && !interp2_analyze(ip))
152 +                return(-1);
153 +                                /* get x & y grid positions */
154 +        fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
155 +
156 +        if (fgi[0] >= NI2DIM)
157 +                fgi[0] = NI2DIM-1;
158 +        else if (fgi[0] < 0)
159 +                fgi[0] = 0;
160 +        else
161 +                ++inbounds;
162 +
163 +        fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
164 +
165 +        if (fgi[1] >= NI2DIM)
166 +                fgi[1] = NI2DIM-1;
167 +        else if (fgi[1] < 0)
168 +                fgi[1] = 0;
169 +        else
170 +                ++inbounds;
171 +
172 +        return(inbounds == 2);
173 + }
174 +
175 + #define setflg(fl,xi,yi)        ((fl)[yi] |= 1<<(xi))
176 +
177 + #define chkflg(fl,xi,yi)        ((fl)[yi]>>(xi) & 1)
178 +
179 + /* private flood function to determine sample influence */
180 + static void
181 + influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
182 +                        int xfi, int yfi)
183 + {
184 +        double  gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
185 +                                        ip->smin[0];
186 +        double  gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
187 +                                        ip->smin[1];
188 +        double  dx = gx - ip->spt[i][0];
189 +        double  dy = gy - ip->spt[i][1];
190 +
191 +        setflg(visited, xfi, yfi);
192 +
193 +        if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
194 +                return;
195 +
196 +        setflg(ip->da[i].infl, xfi, yfi);
197 +
198 +        if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
199 +                influence_flood(ip, i, visited, xfi-1, yfi);
200 +
201 +        if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
202 +                influence_flood(ip, i, visited, xfi, yfi-1);
203 +
204 +        if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
205 +                influence_flood(ip, i, visited, xfi+1, yfi);
206 +
207 +        if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
208 +                influence_flood(ip, i, visited, xfi, yfi+1);
209 + }
210 +
211 + /* private call to compute sample influence maps */
212 + static void
213 + map_influence(INTERP2 *ip)
214 + {
215 +        unsigned short  visited[NI2DIM];
216 +        int             fgi[2];
217 +        int             i, j;
218 +
219 +        for (i = ip->ns; i--; ) {
220 +                for (j = NI2DIM; j--; ) {
221 +                        ip->da[i].infl[j] = 0;
222 +                        visited[j] = 0;
223 +                }
224 +                interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
225 +
226 +                influence_flood(ip, i, visited, fgi[0], fgi[1]);
227 +        }
228 + }
229 +
230   /* Modify smoothing parameter by the given factor */
231   void
232   interp2_smooth(INTERP2 *ip, double sf)
233   {
234 +        float   old_smf = ip->smf;
235 +
236          if ((ip->smf *= sf) < NI2DSMF)
237                  ip->smf = NI2DSMF;
238 +                                        /* need to recompute influence maps? */
239 +        if (ip->da != NULL && (old_smf*.85 > ip->smf) |
240 +                                (ip->smf > old_smf*1.15))
241 +                map_influence(ip);
242   }
243  
244   /* private call-back to sort position index */
# Line 134 | Line 267 | sort_samples(SAMPORD *sord, const INTERP2 *ip, double
267                  sord[i].si = i;
268                  sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
269          }
270 <        qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
270 >        qsort(sord, ip->ns, sizeof(SAMPORD), cmp_spos);
271   }
272  
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 –
273   /* (Re)compute anisotropic basis function interpolant (normally automatic) */
274   int
275   interp2_analyze(INTERP2 *ip)
# Line 170 | Line 290 | interp2_analyze(INTERP2 *ip)
290          ip->smin[0] = ip->smin[1] = FHUGE;
291          ip->smax[0] = ip->smax[1] = -FHUGE;
292          for (i = ip->ns; i--; ) {
293 <                if (ip->spt[i][0] < ip->smin[0])
294 <                        ip->smin[0] = ip->spt[i][0];
295 <                if (ip->spt[i][0] > ip->smax[0])
296 <                        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];
293 >                if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
294 >                if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
295 >                if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
296 >                if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
297          }
298          ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
299                          (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
# Line 185 | Line 301 | interp2_analyze(INTERP2 *ip)
301          if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
302                  return(0);
303                                          /* allocate analysis data */
304 <        ip->da = (struct interp2_samp *)calloc( ip->ns,
305 <                                        sizeof(struct interp2_samp) );
304 >        ip->da = (struct interp2_samp *)malloc(
305 >                                sizeof(struct interp2_samp)*ip->ns );
306          if (ip->da == NULL)
307                  return(0);
308 <                                        /* get temporary arrays */
308 >                                        /* allocate temporary arrays */
309          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
310          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
311          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
# Line 250 | Line 366 | interp2_analyze(INTERP2 *ip)
366                      }
367              }
368          }
369 <        free(sortord);                  /* clean up */
369 >        free(sortord);                  /* release temp arrays */
370          free(rightrndx);
371          free(leftrndx);
372          free(endrndx);
373 <        return(1);
373 >                                        /* map sample influence areas */
374 >        map_influence(ip);
375 >        return(1);                      /* all done */
376   }
377  
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, r2, d2;
265 –        int             ri;
266 –                                /* get relative direction */
267 –        x -= ip->spt[i][0];
268 –        y -= ip->spt[i][1];
269 –        dir = atan2a(y, x);
270 –        dir += 2.*PI*(dir < 0);
271 –                                /* linear radius interpolation */
272 –        rd = dir * (NI2DIR/2./PI);
273 –        ri = (int)rd;
274 –        rd -= (double)ri;
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 –
378   /* Assign full set of normalized weights to interpolate the given position */
379   int
380   interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
381   {
382          double  wnorm;
383          int     fgi[2];
337 –        int     ingrid;
384          int     i;
385  
386          if (wtv == NULL)
387                  return(0);
388                                          /* get flag position */
389 <        if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
389 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
390                  return(0);
391  
392          wnorm = 0;                      /* compute raw weights */
393          for (i = ip->ns; i--; )
394 <            if (chkblk(ip, i, fgi)) {
349 <                wtv[i] = 0;
350 <            } else {
394 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
395                  double  wt = interp2_wti(ip, i, x, y);
396                  wtv[i] = wt;
397                  wnorm += wt;
398 <                if (wt <= 1e-9 && ingrid)
399 <                        setblk(ip, i, fgi);
356 <            }
398 >            } else
399 >                wtv[i] = 0;
400          if (wnorm <= 0)                 /* too far from all our samples! */
401                  return(0);
402          wnorm = 1./wnorm;               /* normalize weights */
# Line 369 | Line 412 | interp2_topsamp(float wt[], int si[], const int n, INT
412   {
413          int     nn = 0;
414          int     fgi[2];
372 –        int     ingrid;
415          double  wnorm;
416          int     i, j;
417  
418          if ((n <= 0) | (wt == NULL) | (si == NULL))
419                  return(0);
420                                          /* get flag position */
421 <        if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
421 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
422                  return(0);
423                                          /* identify top n weights */
424          for (i = ip->ns; i--; )
425 <            if (!chkblk(ip, i, fgi)) {
425 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
426                  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 –                }
427                  for (j = nn; j > 0; j--) {
428                          if (wt[j-1] >= wti)
429                                  break;

Diff Legend

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