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.9 by greg, Tue Feb 12 18:41:39 2013 UTC vs.
Revision 2.10 by greg, Thu Feb 14 19:57:10 2013 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 an array
27 > * of flags indicating where weights are (nearly) zero.
28   ****************************************************************/
29  
30   #include <stdio.h>
# Line 155 | Line 156 | interp2_analyze(INTERP2 *ip)
156   {
157          SAMPORD *sortord;
158          int     *rightrndx, *leftrndx, *endrndx;
159 <        int     bd;
159 >        int     i, bd;
160                                          /* sanity checks */
161 <        if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
161 >        if (ip == NULL)
162                  return(0);
163 <                                        /* need to allocate? */
164 <        if (ip->da == NULL) {
165 <                ip->da = (unsigned short (*)[NI2DIR])malloc(
165 <                                sizeof(unsigned short)*NI2DIR*ip->ns);
166 <                if (ip->da == NULL)
167 <                        return(0);
163 >        if (ip->da != NULL) {           /* free previous data if any */
164 >                free(ip->da);
165 >                ip->da = NULL;
166          }
167 +        if ((ip->ns <= 1) | (ip->dmin <= 0))
168 +                return(0);
169 +                                        /* compute sample domain */
170 +        ip->smin[0] = ip->smin[1] = FHUGE;
171 +        ip->smul[0] = ip->smul[1] = -FHUGE;
172 +        for (i = ip->ns; i--; ) {
173 +                if (ip->spt[i][0] < ip->smin[0])
174 +                        ip->smin[0] = ip->spt[i][0];
175 +                if (ip->spt[i][0] > ip->smul[0])
176 +                        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];
181 +        }
182 +        ip->smul[0] -= ip->smin[0];
183 +        ip->smul[1] -= ip->smin[1];
184 +        ip->grid2 = (ip->smul[0]*ip->smul[0] + ip->smul[1]*ip->smul[1]) *
185 +                        (4./NI2DIM/NI2DIM);
186 +        if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
187 +                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];
192 +                                        /* allocate analysis data */
193 +        ip->da = (struct interp2_samp *)calloc( ip->ns,
194 +                                        sizeof(struct interp2_samp) );
195 +        if (ip->da == NULL)
196 +                return(0);
197                                          /* get temporary arrays */
198          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
199          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
# Line 178 | Line 206 | interp2_analyze(INTERP2 *ip)
206          for (bd = 0; bd < NI2DIR/2; bd++) {
207              const double        ang = 2.*PI/NI2DIR*bd;
208              int                 *sptr;
181 –            int                 i;
209                                          /* create right reverse index */
210              if (bd) {                   /* re-use from previous iteration? */
211                  sptr = rightrndx;
# Line 209 | Line 236 | interp2_analyze(INTERP2 *ip)
236                  const int       ii = sortord[i].si;
237                  int             j;
238                                          /* preload with large radii */
239 <                ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
240 <                            .5*(sortord[ip->ns-1].dm - sortord[0].dm));
239 >                ip->da[ii].dia[bd] =
240 >                ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
241 >                                .5*(sortord[ip->ns-1].dm - sortord[0].dm));
242                  for (j = i; ++j < ip->ns; )     /* nearest above */
243                      if (rightrndx[sortord[j].si] > rightrndx[ii] &&
244                                      leftrndx[sortord[j].si] < leftrndx[ii]) {
245 <                        ip->da[ii][bd] = encode_diameter(ip,
245 >                        ip->da[ii].dia[bd] = encode_diameter(ip,
246                                                  sortord[j].dm - sortord[i].dm);
247                          break;
248                      }
249                  for (j = i; j-- > 0; )          /* nearest below */
250                      if (rightrndx[sortord[j].si] < rightrndx[ii] &&
251                                      leftrndx[sortord[j].si] > leftrndx[ii]) {
252 <                        ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
252 >                        ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
253                                                  sortord[i].dm - sortord[j].dm);
254                          break;
255                      }
# Line 234 | Line 262 | interp2_analyze(INTERP2 *ip)
262          return(1);
263   }
264  
265 < /* private call returns raw weight for a particular sample */
266 < static double
267 < get_wt(const INTERP2 *ip, const int i, double x, double y)
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 <                                /* get relative direction */
273 <        x -= ip->spt[i][0];
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 <        dir = atan2a(y, x);
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][ri] + rd*ip->da[i][(ri+1)%NI2DIR];
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 <        d2 = x*x + y*y;
303 <        if (d2 > 21.*r2)        /* result would be < 1e-9 */
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   }
# Line 268 | Line 316 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
316  
317          if ((wtv == NULL) | (ip == NULL))
318                  return(0);
271 –                                        /* need to compute interpolant? */
272 –        if (ip->da == NULL && !interp2_analyze(ip))
273 –                return(0);
319  
320          wnorm = 0;                      /* compute raw weights */
321          for (i = ip->ns; i--; ) {
322 <                double  wt = get_wt(ip, i, x, y);
322 >                double  wt = interp2_wti(ip, i, x, y);
323                  wtv[i] = wt;
324                  wnorm += wt;
325          }
# Line 297 | Line 342 | interp2_topsamp(float wt[], int si[], const int n, INT
342  
343          if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
344                  return(0);
300 –                                        /* need to compute interpolant? */
301 –        if (ip->da == NULL && !interp2_analyze(ip))
302 –                return(0);
345                                          /* identify top n weights */
346          for (i = ip->ns; i--; ) {
347 <                const double    wti = get_wt(ip, i, x, y);
347 >                const double    wti = interp2_wti(ip, i, x, y);
348                  if (wti <= 1e-9)
349                          continue;
350                  for (j = nn; j > 0; j--) {

Diff Legend

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