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.2 by greg, Sat Feb 9 17:39:21 2013 UTC

# Line 39 | Line 39 | typedef struct {
39          float   dm;             /* distance measure in this direction */
40   } SAMPORD;
41  
42 < /* Allocate a new set of interpolation samples */
42 > /* Allocate a new set of interpolation samples (caller assigns spt[] array) */
43   INTERP2 *
44   interp2_alloc(int nsamps)
45   {
# Line 60 | Line 60 | interp2_alloc(int nsamps)
60          return(nip);
61   }
62  
63 + /* Resize interpolation array (caller must assign any new values) */
64 + INTERP2 *
65 + interp2_realloc(INTERP2 *ip, int nsamps)
66 + {
67 +        if (ip == NULL)
68 +                return(interp2_alloc(nsamps));
69 +        if (nsamps <= 1) {
70 +                interp2_free(ip);
71 +                return(NULL);
72 +        }
73 +        if (nsamps == ip->ns);
74 +                return(ip);
75 +        if (ip->ra != NULL) {   /* will need to recompute distribution */
76 +                free(ip->ra);
77 +                ip->ra = NULL;
78 +        }
79 +        ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
80 +        if (ip == NULL)
81 +                return(NULL);
82 +        ip->ns = nsamps;
83 +        return(ip);
84 + }
85 +
86   /* private call-back to sort position index */
87   static int
88   cmp_spos(const void *p1, const void *p2)
# Line 74 | Line 97 | cmp_spos(const void *p1, const void *p2)
97          return 0;
98   }
99  
100 + /* private routine to order samples in a particular direction */
101 + static void
102 + sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang)
103 + {
104 +        const double    cosd = cos(ang);
105 +        const double    sind = sin(ang);
106 +        int             i;
107 +
108 +        for (i = ip->ns; i--; ) {
109 +                sord[i].si = i;
110 +                sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
111 +        }
112 +        qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
113 + }
114 +
115   /* private routine to encode radius with range checks */
116   static int
117   encode_radius(const INTERP2 *ip, double r)
# Line 87 | Line 125 | encode_radius(const INTERP2 *ip, double r)
125          return(er);
126   }
127  
128 < /* Compute anisotropic Gaussian basis function interpolant */
129 < static int
130 < interp2_compute(INTERP2 *ip)
128 > /* (Re)compute anisotropic basis function interpolant (normally automatic) */
129 > int
130 > interp2_analyze(INTERP2 *ip)
131   {
132          SAMPORD *sortord;
133 <        int     *rightrndx, *leftrndx;
133 >        int     *rightrndx, *leftrndx, *endrndx;
134          int     bd;
135                                          /* sanity checks */
136          if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
# Line 108 | Line 146 | interp2_compute(INTERP2 *ip)
146          sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
147          rightrndx = (int *)malloc(sizeof(int)*ip->ns);
148          leftrndx = (int *)malloc(sizeof(int)*ip->ns);
149 <        if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL))
149 >        endrndx = (int *)malloc(sizeof(int)*ip->ns);
150 >        if ((sortord == NULL) | (rightrndx == NULL) |
151 >                        (leftrndx == NULL) | (endrndx == NULL))
152                  return(0);
153                                          /* run through bidirections */
154          for (bd = 0; bd < NI2DIR/2; bd++) {
155              const double        ang = 2.*PI/NI2DIR*bd;
156 <            double              cosd, sind;
156 >            int                 *sptr;
157              int                 i;
158                                          /* create right reverse index */
159 <            if (bd) {                   /* re-use from prev. iteration? */
160 <                int     *sptr = rightrndx;
159 >            if (bd) {                   /* re-use from previous iteration? */
160 >                sptr = rightrndx;
161                  rightrndx = leftrndx;
162                  leftrndx = sptr;
163 <            } else {                    /* else compute it */
164 <                cosd = cos(ang + (PI/2. - PI/NI2DIR));
165 <                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++)
163 >            } else {                    /* else sort first half-plane */
164 >                sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
165 >                for (i = ip->ns; i--; )
166                      rightrndx[sortord[i].si] = i;
167 +                                        /* & store reverse order for later */
168 +                for (i = ip->ns; i--; )
169 +                    endrndx[sortord[i].si] = ip->ns-1 - i;
170              }
171                                          /* create new left reverse index */
172 <            cosd = cos(ang + (PI/2. + PI/NI2DIR));
173 <            sind = sin(ang + (PI/2. + PI/NI2DIR));
174 <            for (i = 0; i < ip->ns; i++) {
175 <                    sortord[i].si = i;
176 <                    sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
177 <            }
178 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
142 <            for (i = 0; i < ip->ns; i++)
172 >            if (bd == NI2DIR/2 - 1) {   /* use order from first iteration? */
173 >                sptr = leftrndx;
174 >                leftrndx = endrndx;
175 >                endrndx = sptr;
176 >            } else {                    /* else compute new half-plane */
177 >                sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
178 >                for (i = ip->ns; i--; )
179                      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];
180              }
181 <            qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
181 >                                        /* sort grid values in this direction */
182 >            sort_samples(sortord, ip, ang);
183                                          /* find nearest neighbors each side */
184 <            for (i = 0; i < ip->ns; i++) {
184 >            for (i = ip->ns; i--; ) {
185                  const int       rpos = rightrndx[sortord[i].si];
186                  const int       lpos = leftrndx[sortord[i].si];
187                  int             j;
# 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  
# Line 210 | Line 242 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
242          if ((wtv == NULL) | (ip == NULL))
243                  return(0);
244                                          /* need to compute interpolant? */
245 <        if (ip->ra == NULL && !interp2_compute(ip))
245 >        if (ip->ra == NULL && !interp2_analyze(ip))
246                  return(0);
247  
248          wnorm = 0;                      /* compute raw weights */
# Line 244 | Line 276 | interp2_topsamp(float wt[], int si[], const int n, INT
276          if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
277                  return(0);
278                                          /* need to compute interpolant? */
279 <        if (ip->ra == NULL && !interp2_compute(ip))
279 >        if (ip->ra == NULL && !interp2_analyze(ip))
280                  return(0);
281                                          /* identify top n weights */
282          for (i = ip->ns; i--; ) {

Diff Legend

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