ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/common/interp2d.c
(Generate patch)

Comparing src/common/interp2d.c (file contents):
Revision 2.11 by greg, Fri Feb 15 01:26:47 2013 UTC vs.
Revision 2.12 by greg, Fri Feb 15 19:15:16 2013 UTC

# Line 150 | Line 150 | encode_diameter(const INTERP2 *ip, double d)
150          return(ed);
151   }
152  
153 + /* Compute unnormalized weight for a position relative to a sample */
154 + double
155 + interp2_wti(INTERP2 *ip, const int i, double x, double y)
156 + {
157 +        double          dir, rd, r2, d2;
158 +        int             ri;
159 +                                /* get relative direction */
160 +        x -= ip->spt[i][0];
161 +        y -= ip->spt[i][1];
162 +        dir = atan2a(y, x);
163 +        dir += 2.*PI*(dir < 0);
164 +                                /* linear radius interpolation */
165 +        rd = dir * (NI2DIR/2./PI);
166 +        ri = (int)rd;
167 +        rd -= (double)ri;
168 +        rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
169 +        rd = ip->smf * DECODE_DIA(ip, rd);
170 +        r2 = 2.*rd*rd;
171 +        d2 = x*x + y*y;
172 +        if (d2 > 21.*r2)        /* result would be < 1e-9 */
173 +                return(.0);
174 +                                /* Gaussian times harmonic weighting */
175 +        return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
176 + }
177 +
178 + /* private call to get grid flag index */
179 + static int
180 + interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
181 + {
182 +        int     inbounds = 0;
183 +
184 +        if (ip == NULL)         /* paranoia */
185 +                return(-1);
186 +                                /* need to compute interpolant? */
187 +        if (ip->da == NULL && !interp2_analyze(ip))
188 +                return(-1);
189 +                                /* get x & y grid positions */
190 +        fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
191 +
192 +        if (fgi[0] >= NI2DIM)
193 +                fgi[0] = NI2DIM-1;
194 +        else if (fgi[0] < 0)
195 +                fgi[0] = 0;
196 +        else
197 +                ++inbounds;
198 +
199 +        fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
200 +
201 +        if (fgi[1] >= NI2DIM)
202 +                fgi[1] = NI2DIM-1;
203 +        else if (fgi[1] < 0)
204 +                fgi[1] = 0;
205 +        else
206 +                ++inbounds;
207 +
208 +        return(inbounds == 2);
209 + }
210 +
211 + #define setflg(fl,xi,yi)        ((fl)[yi] |= 1<<(xi))
212 +
213 + #define chkflg(fl,xi,yi)        ((fl)[yi]>>(xi) & 1)
214 +
215 + /* private flood function to determine sample influence */
216 + static void
217 + influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
218 +                        int xfi, int yfi)
219 + {
220 +        double  gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
221 +                                        ip->smin[0];
222 +        double  gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
223 +                                        ip->smin[1];
224 +        double  dx = gx - ip->spt[i][0];
225 +        double  dy = gy - ip->spt[i][1];
226 +
227 +        setflg(visited, xfi, yfi);
228 +
229 +        if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
230 +                return;
231 +
232 +        setflg(ip->da[i].infl, xfi, yfi);
233 +
234 +        if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
235 +                influence_flood(ip, i, visited, xfi-1, yfi);
236 +
237 +        if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
238 +                influence_flood(ip, i, visited, xfi, yfi-1);
239 +
240 +        if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
241 +                influence_flood(ip, i, visited, xfi+1, yfi);
242 +
243 +        if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
244 +                influence_flood(ip, i, visited, xfi, yfi+1);
245 + }
246 +
247   /* (Re)compute anisotropic basis function interpolant (normally automatic) */
248   int
249   interp2_analyze(INTERP2 *ip)
250   {
251          SAMPORD *sortord;
252          int     *rightrndx, *leftrndx, *endrndx;
253 <        int     i, bd;
253 >        int     i, j, bd;
254                                          /* sanity checks */
255          if (ip == NULL)
256                  return(0);
# Line 170 | Line 264 | interp2_analyze(INTERP2 *ip)
264          ip->smin[0] = ip->smin[1] = FHUGE;
265          ip->smax[0] = ip->smax[1] = -FHUGE;
266          for (i = ip->ns; i--; ) {
267 <                if (ip->spt[i][0] < ip->smin[0])
268 <                        ip->smin[0] = ip->spt[i][0];
269 <                if (ip->spt[i][0] > ip->smax[0])
270 <                        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];
267 >                if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
268 >                if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
269 >                if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
270 >                if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
271          }
272          ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
273                          (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
# Line 189 | Line 279 | interp2_analyze(INTERP2 *ip)
279                                          sizeof(struct interp2_samp) );
280          if (ip->da == NULL)
281                  return(0);
282 <                                        /* get temporary arrays */
282 >                                        /* allocate temporary arrays */
283          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
284          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
285          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
# Line 229 | Line 319 | interp2_analyze(INTERP2 *ip)
319                                          /* find nearest neighbors each side */
320              for (i = ip->ns; i--; ) {
321                  const int       ii = sortord[i].si;
232 –                int             j;
322                                          /* preload with large radii */
323                  ip->da[ii].dia[bd] =
324                  ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
# Line 250 | Line 339 | interp2_analyze(INTERP2 *ip)
339                      }
340              }
341          }
342 <        free(sortord);                  /* clean up */
342 >        free(sortord);                  /* release temp arrays */
343          free(rightrndx);
344          free(leftrndx);
345          free(endrndx);
346 <        return(1);
347 < }
346 >                                        /* fill influence maps */
347 >        for (i = ip->ns; i--; ) {
348 >                unsigned short  visited[NI2DIM];
349 >                int             fgi[2];
350  
351 < /* Compute unnormalized weight for a position relative to a sample */
352 < double
353 < 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;
351 >                for (j = NI2DIM; j--; ) visited[j] = 0;
352 >                interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
353 >                influence_flood(ip, i, visited, fgi[0], fgi[1]);
354          }
355 <        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);
355 >        return(1);                      /* all done */
356   }
357  
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 –
358   /* Assign full set of normalized weights to interpolate the given position */
359   int
360   interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
361   {
362          double  wnorm;
363          int     fgi[2];
337 –        int     ingrid;
364          int     i;
365  
366          if (wtv == NULL)
367                  return(0);
368                                          /* get flag position */
369 <        if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
369 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
370                  return(0);
371  
372          wnorm = 0;                      /* compute raw weights */
373          for (i = ip->ns; i--; )
374 <            if (chkblk(ip, i, fgi)) {
349 <                wtv[i] = 0;
350 <            } else {
374 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
375                  double  wt = interp2_wti(ip, i, x, y);
376                  wtv[i] = wt;
377                  wnorm += wt;
378 <                if (wt <= 1e-9 && ingrid)
379 <                        setblk(ip, i, fgi);
356 <            }
378 >            } else
379 >                wtv[i] = 0;
380          if (wnorm <= 0)                 /* too far from all our samples! */
381                  return(0);
382          wnorm = 1./wnorm;               /* normalize weights */
# Line 369 | Line 392 | interp2_topsamp(float wt[], int si[], const int n, INT
392   {
393          int     nn = 0;
394          int     fgi[2];
372 –        int     ingrid;
395          double  wnorm;
396          int     i, j;
397  
398          if ((n <= 0) | (wt == NULL) | (si == NULL))
399                  return(0);
400                                          /* get flag position */
401 <        if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
401 >        if (interp2_flagpos(fgi, ip, x, y) < 0)
402                  return(0);
403                                          /* identify top n weights */
404          for (i = ip->ns; i--; )
405 <            if (!chkblk(ip, i, fgi)) {
405 >            if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
406                  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 –                }
407                  for (j = nn; j > 0; j--) {
408                          if (wt[j-1] >= wti)
409                                  break;

Diff Legend

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