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.1 by greg, Sat Feb 9 00:55:40 2013 UTC vs.
Revision 2.3 by greg, Mon Feb 11 20:01:15 2013 UTC

# Line 18 | Line 18 | static const char RCSid[] = "$Id$";
18   * calculation, we sort the data into 3 half-planes and
19   * perform simple tests to see which neighbor is closest in
20   * a each direction.  Once we have our approximate neighborhood
21 < * for a sample, we can use it in a Gaussian weighting scheme
22 < * with anisotropic surround.  This gives us a fairly smooth
23 < * interpolation however the sample points may be initially
24 < * distributed.  Evaluation is accelerated by use of a fast
25 < * approximation to the atan2(y,x) function.
21 > * for a sample, we can use it in a modified Gaussian weighting
22 > * scheme with anisotropic surround.  Harmonic weighting is added
23 > * to reduce the influence of distant neighbors.  This yields a
24 > * smooth interpolation regardless of how the sample points are
25 > * initiallydistributed.  Evaluation is accelerated by use of a
26 > * fast approximation to the atan2(y,x) function.
27   **************************************************************/
28  
29   #include <stdio.h>
# Line 39 | Line 40 | typedef struct {
40          float   dm;             /* distance measure in this direction */
41   } SAMPORD;
42  
43 < /* Allocate a new set of interpolation samples */
43 > /* Allocate a new set of interpolation samples (caller assigns spt[] array) */
44   INTERP2 *
45   interp2_alloc(int nsamps)
46   {
# Line 60 | Line 61 | interp2_alloc(int nsamps)
61          return(nip);
62   }
63  
64 + /* Resize interpolation array (caller must assign any new values) */
65 + INTERP2 *
66 + interp2_realloc(INTERP2 *ip, int nsamps)
67 + {
68 +        if (ip == NULL)
69 +                return(interp2_alloc(nsamps));
70 +        if (nsamps <= 1) {
71 +                interp2_free(ip);
72 +                return(NULL);
73 +        }
74 +        if (nsamps == ip->ns);
75 +                return(ip);
76 +        if (ip->ra != NULL) {   /* will need to recompute distribution */
77 +                free(ip->ra);
78 +                ip->ra = NULL;
79 +        }
80 +        ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
81 +        if (ip == NULL)
82 +                return(NULL);
83 +        ip->ns = nsamps;
84 +        return(ip);
85 + }
86 +
87   /* private call-back to sort position index */
88   static int
89   cmp_spos(const void *p1, const void *p2)
# Line 74 | Line 98 | cmp_spos(const void *p1, const void *p2)
98          return 0;
99   }
100  
101 + /* private routine to order samples in a particular direction */
102 + static void
103 + sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang)
104 + {
105 +        const double    cosd = cos(ang);
106 +        const double    sind = sin(ang);
107 +        int             i;
108 +
109 +        for (i = ip->ns; i--; ) {
110 +                sord[i].si = i;
111 +                sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
112 +        }
113 +        qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
114 + }
115 +
116   /* private routine to encode radius with range checks */
117   static int
118   encode_radius(const INTERP2 *ip, double r)
# Line 87 | Line 126 | encode_radius(const INTERP2 *ip, double r)
126          return(er);
127   }
128  
129 < /* Compute anisotropic Gaussian basis function interpolant */
130 < static int
131 < interp2_compute(INTERP2 *ip)
129 > /* (Re)compute anisotropic basis function interpolant (normally automatic) */
130 > int
131 > interp2_analyze(INTERP2 *ip)
132   {
133          SAMPORD *sortord;
134 <        int     *rightrndx, *leftrndx;
134 >        int     *rightrndx, *leftrndx, *endrndx;
135          int     bd;
136                                          /* sanity checks */
137          if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
# Line 108 | Line 147 | interp2_compute(INTERP2 *ip)
147          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
148          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
149          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
150 <        if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL))
150 >        endrndx = (int *)malloc(sizeof(int)*ip->ns);
151 >        if ((sortord == NULL) | (rightrndx == NULL) |
152 >                        (leftrndx == NULL) | (endrndx == NULL))
153                  return(0);
154                                          /* run through bidirections */
155          for (bd = 0; bd < NI2DIR/2; bd++) {
156              const double        ang = 2.*PI/NI2DIR*bd;
157 <            double              cosd, sind;
157 >            int                 *sptr;
158              int                 i;
159                                          /* create right reverse index */
160 <            if (bd) {                   /* re-use from prev. iteration? */
161 <                int     *sptr = rightrndx;
160 >            if (bd) {                   /* re-use from previous iteration? */
161 >                sptr = rightrndx;
162                  rightrndx = leftrndx;
163                  leftrndx = sptr;
164 <            } else {                    /* else compute it */
165 <                cosd = cos(ang + (PI/2. - PI/NI2DIR));
166 <                sind = sin(ang + (PI/2. - PI/NI2DIR));
126 <                for (i = 0; i < ip->ns; i++) {
127 <                    sortord[i].si = i;
128 <                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
129 <                }
130 <                qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
131 <                for (i = 0; i < ip->ns; i++)
164 >            } else {                    /* else sort first half-plane */
165 >                sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
166 >                for (i = ip->ns; i--; )
167                      rightrndx[sortord[i].si] = i;
168 +                                        /* & store reverse order for later */
169 +                for (i = ip->ns; i--; )
170 +                    endrndx[sortord[i].si] = ip->ns-1 - i;
171              }
172                                          /* create new left reverse index */
173 <            cosd = cos(ang + (PI/2. + PI/NI2DIR));
174 <            sind = sin(ang + (PI/2. + PI/NI2DIR));
175 <            for (i = 0; i < ip->ns; i++) {
176 <                    sortord[i].si = i;
177 <                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
178 <            }
179 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
142 <            for (i = 0; i < ip->ns; i++)
173 >            if (bd == NI2DIR/2 - 1) {   /* use order from first iteration? */
174 >                sptr = leftrndx;
175 >                leftrndx = endrndx;
176 >                endrndx = sptr;
177 >            } else {                    /* else compute new half-plane */
178 >                sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
179 >                for (i = ip->ns; i--; )
180                      leftrndx[sortord[i].si] = i;
144 –                                        /* sort grid values in this direction */
145 –            cosd = cos(ang);
146 –            sind = sin(ang);
147 –            for (i = 0; i < ip->ns; i++) {
148 –                    sortord[i].si = i;
149 –                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
181              }
182 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
182 >                                        /* sort grid values in this direction */
183 >            sort_samples(sortord, ip, ang);
184                                          /* find nearest neighbors each side */
185 <            for (i = 0; i < ip->ns; i++) {
186 <                const int       rpos = rightrndx[sortord[i].si];
155 <                const int       lpos = leftrndx[sortord[i].si];
185 >            for (i = ip->ns; i--; ) {
186 >                const int       ii = sortord[i].si;
187                  int             j;
188 <                                        /* preload with large radius */
189 <                ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
188 >                                        /* preload with large radii */
189 >                ip->ra[ii][bd] = ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
190                              .25*(sortord[ip->ns-1].dm - sortord[0].dm));
191                  for (j = i; ++j < ip->ns; )     /* nearest above */
192 <                    if (rightrndx[sortord[j].si] > rpos &&
193 <                                    leftrndx[sortord[j].si] < lpos) {
194 <                        ip->ra[i][bd] = encode_radius(ip,
192 >                    if (rightrndx[sortord[j].si] > rightrndx[ii] &&
193 >                                    leftrndx[sortord[j].si] < leftrndx[ii]) {
194 >                        ip->ra[ii][bd] = encode_radius(ip,
195                                          .5*(sortord[j].dm - sortord[i].dm));
196                          break;
197                      }
198                  for (j = i; j-- > 0; )          /* nearest below */
199 <                    if (rightrndx[sortord[j].si] < rpos &&
200 <                                    leftrndx[sortord[j].si] > lpos) {
201 <                        ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
199 >                    if (rightrndx[sortord[j].si] < rightrndx[ii] &&
200 >                                    leftrndx[sortord[j].si] > leftrndx[ii]) {
201 >                        ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
202                                          .5*(sortord[i].dm - sortord[j].dm));
203                          break;
204                      }
# Line 176 | Line 207 | interp2_compute(INTERP2 *ip)
207          free(sortord);                  /* clean up */
208          free(rightrndx);
209          free(leftrndx);
210 +        free(endrndx);
211          return(1);
212   }
213  
214 < /* private call returns log of raw weight for a particular sample */
214 > /* private call returns raw weight for a particular sample */
215   static double
216 < get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
216 > get_wt(const INTERP2 *ip, const int i, double x, double y)
217   {
218 <        double  dir, rd;
218 >        double  dir, rd, d2;
219          int     ri;
220                                  /* get relative direction */
221          x -= ip->spt[i][0];
# Line 196 | Line 228 | get_ln_wt(const INTERP2 *ip, const int i, double x, do
228          rd -= (double)ri;
229          rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
230          rd = ip->smf * DECODE_RAD(ip, rd);
231 <                                /* return log of Gaussian weight */
232 <        return( (x*x + y*y) / (-2.*rd*rd) );
231 >        d2 = x*x + y*y;
232 >                                /* Gaussian times harmonic weighting */
233 >        return( exp(d2/(-2.*rd*rd)) * ip->rmin/(ip->rmin + sqrt(d2)) );
234   }
235  
236   /* Assign full set of normalized weights to interpolate the given position */
# Line 210 | Line 243 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
243          if ((wtv == NULL) | (ip == NULL))
244                  return(0);
245                                          /* need to compute interpolant? */
246 <        if (ip->ra == NULL && !interp2_compute(ip))
246 >        if (ip->ra == NULL && !interp2_analyze(ip))
247                  return(0);
248  
249          wnorm = 0;                      /* compute raw weights */
250          for (i = ip->ns; i--; ) {
251 <                double  wt = get_ln_wt(ip, i, x, y);
219 <                if (wt < -21.) {
220 <                        wtv[i] = 0;     /* ignore weights < 1e-9 */
221 <                        continue;
222 <                }
223 <                wt = exp(wt);           /* Gaussian weight */
251 >                double  wt = get_wt(ip, i, x, y);
252                  wtv[i] = wt;
253                  wnorm += wt;
254          }
# Line 244 | Line 272 | interp2_topsamp(float wt[], int si[], const int n, INT
272          if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
273                  return(0);
274                                          /* need to compute interpolant? */
275 <        if (ip->ra == NULL && !interp2_compute(ip))
275 >        if (ip->ra == NULL && !interp2_analyze(ip))
276                  return(0);
277                                          /* identify top n weights */
278          for (i = ip->ns; i--; ) {
279 <                const double    lnwt = get_ln_wt(ip, i, x, y);
279 >                const double    wti = get_wt(ip, i, x, y);
280                  for (j = nn; j > 0; j--) {
281 <                        if (wt[j-1] >= lnwt)
281 >                        if (wt[j-1] >= wti)
282                                  break;
283                          if (j < n) {
284                                  wt[j] = wt[j-1];
# Line 258 | Line 286 | interp2_topsamp(float wt[], int si[], const int n, INT
286                          }
287                  }
288                  if (j < n) {            /* add/insert sample */
289 <                        wt[j] = lnwt;
289 >                        wt[j] = wti;
290                          si[j] = i;
291                          nn += (nn < n);
292                  }
293          }
294 <        wnorm = 0;                      /* exponentiate and normalize */
295 <        for (j = nn; j--; ) {
296 <                double  dwt = exp(wt[j]);
269 <                wt[j] = dwt;
270 <                wnorm += dwt;
271 <        }
294 >        wnorm = 0;                      /* normalize sample weights */
295 >        for (j = nn; j--; )
296 >                wnorm += wt[j];
297          if (wnorm <= 0)
298                  return(0);
299          wnorm = 1./wnorm;

Diff Legend

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