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

Comparing ray/src/cv/bsdfinterp.c (file contents):
Revision 2.1 by greg, Fri Oct 19 04:14:29 2012 UTC vs.
Revision 2.12 by greg, Thu Sep 26 17:05:00 2013 UTC

# Line 13 | Line 13 | static const char RCSid[] = "$Id$";
13   #include <string.h>
14   #include <math.h>
15   #include "bsdfrep.h"
16 –                                /* migration edges drawn in raster fashion */
17 – MIGRATION               *mig_grid[GRIDRES][GRIDRES];
16  
19 – #ifdef DEBUG
20 – #include "random.h"
21 – #include "bmpfile.h"
22 – /* Hash pointer to byte value (must return 0 for NULL) */
23 – static int
24 – byte_hash(const void *p)
25 – {
26 –        size_t  h = (size_t)p;
27 –        h ^= (size_t)p >> 8;
28 –        h ^= (size_t)p >> 16;
29 –        h ^= (size_t)p >> 24;
30 –        return(h & 0xff);
31 – }
32 – /* Write out BMP image showing edges */
33 – static void
34 – write_edge_image(const char *fname)
35 – {
36 –        BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
37 –        BMPWriter       *wtr;
38 –        int             i, j;
39 –
40 –        fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
41 –        hdr->compr = BI_RLE8;
42 –        for (i = 256; --i; ) {                  /* assign random color map */
43 –                hdr->palette[i].r = random() & 0xff;
44 –                hdr->palette[i].g = random() & 0xff;
45 –                hdr->palette[i].b = random() & 0xff;
46 –                                                /* reject dark colors */
47 –                i += (hdr->palette[i].r + hdr->palette[i].g +
48 –                                                hdr->palette[i].b < 128);
49 –        }
50 –        hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
51 –                                                /* open output */
52 –        wtr = BMPopenOutputFile(fname, hdr);
53 –        if (wtr == NULL) {
54 –                free(hdr);
55 –                return;
56 –        }
57 –        for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
58 –                for (j = 0; j < GRIDRES; j++)
59 –                        wtr->scanline[j] = byte_hash(mig_grid[i][j]);
60 –                if (BMPwriteScanline(wtr) != BIR_OK)
61 –                        break;
62 –        }
63 –        BMPcloseOutput(wtr);                    /* close & clean up */
64 – }
65 – #endif
66 –
67 – /* Draw edge list into mig_grid array */
68 – void
69 – draw_edges(void)
70 – {
71 –        int             nnull = 0, ntot = 0;
72 –        MIGRATION       *ej;
73 –        int             p0[2], p1[2];
74 –
75 –        memset(mig_grid, 0, sizeof(mig_grid));
76 –        for (ej = mig_list; ej != NULL; ej = ej->next) {
77 –                ++ntot;
78 –                pos_from_vec(p0, ej->rbfv[0]->invec);
79 –                pos_from_vec(p1, ej->rbfv[1]->invec);
80 –                if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
81 –                        ++nnull;
82 –                        mig_grid[p0[0]][p0[1]] = ej;
83 –                        continue;
84 –                }
85 –                if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
86 –                        const int       xstep = 2*(p1[0] > p0[0]) - 1;
87 –                        const double    ystep = (double)((p1[1]-p0[1])*xstep) /
88 –                                                        (double)(p1[0]-p0[0]);
89 –                        int             x;
90 –                        double          y;
91 –                        for (x = p0[0], y = p0[1]+.5; x != p1[0];
92 –                                                x += xstep, y += ystep)
93 –                                mig_grid[x][(int)y] = ej;
94 –                        mig_grid[x][(int)y] = ej;
95 –                } else {
96 –                        const int       ystep = 2*(p1[1] > p0[1]) - 1;
97 –                        const double    xstep = (double)((p1[0]-p0[0])*ystep) /
98 –                                                        (double)(p1[1]-p0[1]);
99 –                        int             y;
100 –                        double          x;
101 –                        for (y = p0[1], x = p0[0]+.5; y != p1[1];
102 –                                                y += ystep, x += xstep)
103 –                                mig_grid[(int)x][y] = ej;
104 –                        mig_grid[(int)x][y] = ej;
105 –                }
106 –        }
107 –        if (nnull)
108 –                fprintf(stderr, "Warning: %d of %d edges are null\n",
109 –                                nnull, ntot);
110 – #ifdef DEBUG
111 –        write_edge_image("bsdf_edges.bmp");
112 – #endif
113 – }
114 –
115 – /* Identify enclosing triangle for this position (flood fill raster check) */
116 – static int
117 – identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
118 –                        int px, int py)
119 – {
120 –        const int       btest = 1<<(py&07);
121 –
122 –        if (vmap[px][py>>3] & btest)            /* already visited here? */
123 –                return(1);
124 –                                                /* else mark it */
125 –        vmap[px][py>>3] |= btest;
126 –
127 –        if (mig_grid[px][py] != NULL) {         /* are we on an edge? */
128 –                int     i;
129 –                for (i = 0; i < 3; i++) {
130 –                        if (miga[i] == mig_grid[px][py])
131 –                                return(1);
132 –                        if (miga[i] != NULL)
133 –                                continue;
134 –                        miga[i] = mig_grid[px][py];
135 –                        return(1);
136 –                }
137 –                return(0);                      /* outside triangle! */
138 –        }
139 –                                                /* check neighbors (flood) */
140 –        if (px > 0 && !identify_tri(miga, vmap, px-1, py))
141 –                return(0);
142 –        if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
143 –                return(0);
144 –        if (py > 0 && !identify_tri(miga, vmap, px, py-1))
145 –                return(0);
146 –        if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
147 –                return(0);
148 –        return(1);                              /* this neighborhood done */
149 – }
150 –
17   /* Insert vertex in ordered list */
18   static void
19   insert_vert(RBFNODE **vlist, RBFNODE *v)
# Line 185 | Line 51 | order_triangle(MIGRATION *miga[3])
51                  insert_vert(vert, miga[i]->rbfv[1]);
52          }
53                                                  /* should be just 3 vertices */
54 <        if ((vert[3] == NULL) | (vert[4] != NULL))
54 >        if ((vert[2] == NULL) | (vert[3] != NULL))
55                  return(0);
56                                                  /* identify edge 0 */
57          for (i = 3; i--; )
# Line 219 | Line 85 | order_triangle(MIGRATION *miga[3])
85          return(1);
86   }
87  
88 + /* Determine if we are close enough to an edge */
89 + static int
90 + on_edge(const MIGRATION *ej, const FVECT ivec)
91 + {
92 +        double  cos_a, cos_b, cos_c, cos_aplusb;
93 +                                        /* use triangle inequality */
94 +        cos_a = DOT(ej->rbfv[0]->invec, ivec);
95 +        if (cos_a <= 0)
96 +                return(0);
97 +
98 +        cos_b = DOT(ej->rbfv[1]->invec, ivec);
99 +        if (cos_b <= 0)
100 +                return(0);
101 +
102 +        cos_aplusb = cos_a*cos_b - sqrt((1.-cos_a*cos_a)*(1.-cos_b*cos_b));
103 +        if (cos_aplusb <= 0)
104 +                return(0);
105 +
106 +        cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec);
107 +
108 +        return(cos_c - cos_aplusb < .001);
109 + }
110 +
111 + /* Determine if we are inside the given triangle */
112 + static int
113 + in_tri(const RBFNODE *v1, const RBFNODE *v2, const RBFNODE *v3, const FVECT p)
114 + {
115 +        FVECT   vc;
116 +        int     sgn1, sgn2, sgn3;
117 +                                        /* signed volume test */
118 +        VCROSS(vc, v1->invec, v2->invec);
119 +        sgn1 = (DOT(p, vc) > 0);
120 +        VCROSS(vc, v2->invec, v3->invec);
121 +        sgn2 = (DOT(p, vc) > 0);
122 +        if (sgn1 != sgn2)
123 +                return(0);
124 +        VCROSS(vc, v3->invec, v1->invec);
125 +        sgn3 = (DOT(p, vc) > 0);
126 +        return(sgn2 == sgn3);
127 + }
128 +
129 + /* Test (and set) bitmap for edge */
130 + static int
131 + check_edge(unsigned char *emap, int nedges, const MIGRATION *mig, int mark)
132 + {
133 +        int     ejndx, bit2check;
134 +
135 +        if (mig->rbfv[0]->ord > mig->rbfv[1]->ord)
136 +                ejndx = mig->rbfv[1]->ord + (nedges-1)*mig->rbfv[0]->ord;
137 +        else
138 +                ejndx = mig->rbfv[0]->ord + (nedges-1)*mig->rbfv[1]->ord;
139 +
140 +        bit2check = 1<<(ejndx&07);
141 +
142 +        if (emap[ejndx>>3] & bit2check)
143 +                return(0);
144 +        if (mark)
145 +                emap[ejndx>>3] |= bit2check;
146 +        return(1);
147 + }
148 +
149 + /* Compute intersection with the given position over remaining mesh */
150 + static int
151 + in_mesh(MIGRATION *miga[3], unsigned char *emap, int nedges,
152 +                        const FVECT ivec, MIGRATION *mig)
153 + {
154 +        RBFNODE         *tv[2];
155 +        MIGRATION       *sej[2], *dej[2];
156 +        int             i;
157 +                                                /* check visitation record */
158 +        if (!check_edge(emap, nedges, mig, 1))
159 +                return(0);
160 +        if (on_edge(mig, ivec)) {
161 +                miga[0] = mig;                  /* close enough to edge */
162 +                return(1);
163 +        }
164 +        if (!get_triangles(tv, mig))            /* do triangles either side? */
165 +                return(0);
166 +        for (i = 2; i--; ) {                    /* identify edges to check */
167 +                MIGRATION       *ej;
168 +                sej[i] = dej[i] = NULL;
169 +                if (tv[i] == NULL)
170 +                        continue;
171 +                for (ej = tv[i]->ejl; ej != NULL; ej = nextedge(tv[i],ej)) {
172 +                        RBFNODE *rbfop = opp_rbf(tv[i],ej);
173 +                        if (rbfop == mig->rbfv[0]) {
174 +                                if (check_edge(emap, nedges, ej, 0))
175 +                                        sej[i] = ej;
176 +                        } else if (rbfop == mig->rbfv[1]) {
177 +                                if (check_edge(emap, nedges, ej, 0))
178 +                                        dej[i] = ej;
179 +                        }
180 +                }
181 +        }
182 +        for (i = 2; i--; ) {                    /* check triangles just once */
183 +                if (sej[i] != NULL && in_mesh(miga, emap, nedges, ivec, sej[i]))
184 +                        return(1);
185 +                if (dej[i] != NULL && in_mesh(miga, emap, nedges, ivec, dej[i]))
186 +                        return(1);
187 +                if ((sej[i] == NULL) | (dej[i] == NULL))
188 +                        continue;
189 +                if (in_tri(mig->rbfv[0], mig->rbfv[1], tv[i], ivec)) {
190 +                        miga[0] = mig;
191 +                        miga[1] = sej[i];
192 +                        miga[2] = dej[i];
193 +                        return(1);
194 +                }
195 +        }
196 +        return(0);                              /* not near this edge */
197 + }
198 +
199   /* Find edge(s) for interpolating the given vector, applying symmetry */
200   int
201   get_interp(MIGRATION *miga[3], FVECT invec)
202   {
203          miga[0] = miga[1] = miga[2] = NULL;
204          if (single_plane_incident) {            /* isotropic BSDF? */
205 <                RBFNODE *rbf;                   /* find edge we're on */
206 <                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
207 <                        if (input_orient*rbf->invec[2] < input_orient*invec[2])
208 <                                break;
209 <                        if (rbf->next != NULL &&
233 <                                        input_orient*rbf->next->invec[2] <
205 >            RBFNODE     *rbf;                   /* find edge we're on */
206 >            for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
207 >                if (input_orient*rbf->invec[2] < input_orient*invec[2])
208 >                        break;
209 >                if (rbf->next != NULL && input_orient*rbf->next->invec[2] <
210                                                          input_orient*invec[2]) {
211 <                                for (miga[0] = rbf->ejl; miga[0] != NULL;
212 <                                                miga[0] = nextedge(rbf,miga[0]))
213 <                                        if (opp_rbf(rbf,miga[0]) == rbf->next)
214 <                                                return(0);
215 <                                break;
211 >                    for (miga[0] = rbf->ejl; miga[0] != NULL;
212 >                                        miga[0] = nextedge(rbf,miga[0]))
213 >                        if (opp_rbf(rbf,miga[0]) == rbf->next) {
214 >                                double  nf = 1. - rbf->invec[2]*rbf->invec[2];
215 >                                if (nf > FTINY) {       /* rotate to match */
216 >                                        nf = sqrt((1.-invec[2]*invec[2])/nf);
217 >                                        invec[0] = nf*rbf->invec[0];
218 >                                        invec[1] = nf*rbf->invec[1];
219 >                                }
220 >                                return(0);
221                          }
222 +                    break;
223                  }
224 <                return(-1);                     /* outside range! */
224 >            }
225 >            return(-1);                         /* outside range! */
226          }
227          {                                       /* else use triangle mesh */
228 <                const int       sym = use_symmetry(invec);
229 <                unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
230 <                int             pstart[2];
231 <                RBFNODE         *vother;
232 <                MIGRATION       *ej;
233 <                int             i;
234 <
235 <                pos_from_vec(pstart, invec);
236 <                memset(floodmap, 0, sizeof(floodmap));
237 <                                                /* call flooding function */
238 <                if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
239 <                        return(-1);             /* outside mesh */
240 <                if ((miga[0] == NULL) | (miga[2] == NULL))
241 <                        return(-1);             /* should never happen */
242 <                if (miga[1] == NULL)
260 <                        return(sym);            /* on edge */
261 <                                                /* verify triangle */
262 <                if (!order_triangle(miga)) {
228 >                int             sym = use_symmetry(invec);
229 >                int             nedges = 0;
230 >                MIGRATION       *mep;
231 >                unsigned char   *emap;
232 >                                                /* clear visitation map */
233 >                for (mep = mig_list; mep != NULL; mep = mep->next)
234 >                        ++nedges;
235 >                emap = (unsigned char *)calloc((nedges*(nedges-1) + 7)>>3, 1);
236 >                if (emap == NULL) {
237 >                        fprintf(stderr, "%s: Out of memory in get_interp()\n",
238 >                                        progname);
239 >                        exit(1);
240 >                }
241 >                                                /* identify intersection  */
242 >                if (!in_mesh(miga, emap, nedges, invec, mig_list)) {
243   #ifdef DEBUG
244 <                        fputs("Munged triangle in get_interp()\n", stderr);
244 >                        fprintf(stderr,
245 >                        "Incident angle (%.1f,%.1f) deg. outside mesh\n",
246 >                                        get_theta180(invec), get_phi360(invec));
247   #endif
248 <                        vother = NULL;          /* find triangle from edge */
249 <                        for (i = 3; i--; ) {
250 <                            RBFNODE     *tpair[2];
269 <                            if (get_triangles(tpair, miga[i]) &&
270 <                                        (vother = tpair[ is_rev_tri(
271 <                                                        miga[i]->rbfv[0]->invec,
272 <                                                        miga[i]->rbfv[1]->invec,
273 <                                                        invec) ]) != NULL)
274 <                                        break;
275 <                        }
276 <                        if (vother == NULL) {   /* couldn't find 3rd vertex */
248 >                        sym = -1;               /* outside mesh */
249 >                } else if (miga[1] != NULL &&
250 >                                (miga[2] == NULL || !order_triangle(miga))) {
251   #ifdef DEBUG
252 <                                fputs("No triangle in get_interp()\n", stderr);
252 >                        fputs("Munged triangle in get_interp()\n", stderr);
253   #endif
254 <                                return(-1);
281 <                        }
282 <                                                /* reassign other two edges */
283 <                        for (ej = vother->ejl; ej != NULL;
284 <                                                ej = nextedge(vother,ej)) {
285 <                                RBFNODE *vorig = opp_rbf(vother,ej);
286 <                                if (vorig == miga[i]->rbfv[0])
287 <                                        miga[(i+1)%3] = ej;
288 <                                else if (vorig == miga[i]->rbfv[1])
289 <                                        miga[(i+2)%3] = ej;
290 <                        }
291 <                        if (!order_triangle(miga)) {
292 < #ifdef DEBUG
293 <                                fputs("Bad triangle in get_interp()\n", stderr);
294 < #endif
295 <                                return(-1);
296 <                        }
254 >                        sym = -1;
255                  }
256 +                free(emap);
257                  return(sym);                    /* return in standard order */
258          }
259   }
# Line 307 | Line 266 | e_advect_rbf(const MIGRATION *mig, const FVECT invec)
266          int             n, i, j;
267          double          t, full_dist;
268                                                  /* get relative position */
269 <        t = acos(DOT(invec, mig->rbfv[0]->invec));
270 <        if (t < M_PI/GRIDRES) {                 /* near first DSF */
269 >        t = Acos(DOT(invec, mig->rbfv[0]->invec));
270 >        if (t < M_PI/grid_res) {                /* near first DSF */
271                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
272                  rbf = (RBFNODE *)malloc(n);
273                  if (rbf == NULL)
274                          goto memerr;
275                  memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
276 +                rbf->next = NULL; rbf->ejl = NULL;
277                  return(rbf);
278          }
279          full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
280 <        if (t > full_dist-M_PI/GRIDRES) {       /* near second DSF */
280 >        if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
281                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
282                  rbf = (RBFNODE *)malloc(n);
283                  if (rbf == NULL)
284                          goto memerr;
285                  memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
286 +                rbf->next = NULL; rbf->ejl = NULL;
287                  return(rbf);
288          }
289          t /= full_dist;
290          n = 0;                                  /* count migrating particles */
291          for (i = 0; i < mtx_nrows(mig); i++)
292              for (j = 0; j < mtx_ncols(mig); j++)
293 <                n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
293 >                n += (mtx_coef(mig,i,j) > FTINY);
294   #ifdef DEBUG
295          fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
296                          mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
# Line 350 | Line 311 | e_advect_rbf(const MIGRATION *mig, const FVECT invec)
311              float               mv;
312              ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
313              for (j = 0; j < mtx_ncols(mig); j++)
314 <                if ((mv = mig->mtx[mtx_ndx(mig,i,j)]) > FTINY) {
314 >                if ((mv = mtx_coef(mig,i,j)) > FTINY) {
315                          const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
316                          double          rad1 = R2ANG(rbf1j->crad);
317                          FVECT           v;
# Line 376 | Line 337 | memerr:
337  
338   /* Partially advect between recorded incident angles and allocate new RBF */
339   RBFNODE *
340 < advect_rbf(const FVECT invec)
340 > advect_rbf(const FVECT invec, int lobe_lim)
341   {
342 +        double          cthresh = FTINY;
343          FVECT           sivec;
344          MIGRATION       *miga[3];
345          RBFNODE         *rbf;
# Line 393 | Line 355 | advect_rbf(const FVECT invec)
355                  return(NULL);
356          if (miga[1] == NULL) {                  /* advect along edge? */
357                  rbf = e_advect_rbf(miga[0], sivec);
358 <                rev_rbf_symmetry(rbf, sym);
358 >                if (single_plane_incident)
359 >                        rotate_rbf(rbf, invec);
360 >                else
361 >                        rev_rbf_symmetry(rbf, sym);
362                  return(rbf);
363          }
364   #ifdef DEBUG
# Line 415 | Line 380 | advect_rbf(const FVECT invec)
380          geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
381                          s, GEOD_REL);
382          t = acos(DOT(v1,sivec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
383 + tryagain:
384          n = 0;                                  /* count migrating particles */
385          for (i = 0; i < mtx_nrows(miga[0]); i++)
386              for (j = 0; j < mtx_ncols(miga[0]); j++)
387 <                for (k = (miga[0]->mtx[mtx_ndx(miga[0],i,j)] > FTINY) *
387 >                for (k = (mtx_coef(miga[0],i,j) > cthresh) *
388                                          mtx_ncols(miga[2]); k--; )
389 <                        n += (miga[2]->mtx[mtx_ndx(miga[2],i,k)] > FTINY &&
390 <                                miga[1]->mtx[mtx_ndx(miga[1],j,k)] > FTINY);
389 >                        n += (mtx_coef(miga[2],i,k) > cthresh ||
390 >                                mtx_coef(miga[1],j,k) > cthresh);
391 >        if ((lobe_lim > 0) & (n > lobe_lim)) {
392 >                cthresh = cthresh*2. + 10.*FTINY;
393 >                goto tryagain;
394 >        }
395   #ifdef DEBUG
396          fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
397                          miga[0]->rbfv[0]->nrbf, miga[0]->rbfv[1]->nrbf,
# Line 446 | Line 416 | advect_rbf(const FVECT invec)
416              const double        rad0i = R2ANG(rbf0i->crad);
417              ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
418              for (j = 0; j < mtx_ncols(miga[0]); j++) {
419 <                const float     ma = miga[0]->mtx[mtx_ndx(miga[0],i,j)];
419 >                const float     ma = mtx_coef(miga[0],i,j);
420                  const RBFVAL    *rbf1j;
421                  double          rad1j, srad2;
422 <                if (ma <= FTINY)
422 >                if (ma <= cthresh)
423                          continue;
424                  rbf1j = &miga[0]->rbfv[1]->rbfa[j];
425                  rad1j = R2ANG(rbf1j->crad);
# Line 457 | Line 427 | advect_rbf(const FVECT invec)
427                  ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
428                  geodesic(v1, v0, v1, s, GEOD_REL);
429                  for (k = 0; k < mtx_ncols(miga[2]); k++) {
430 <                    float               mb = miga[1]->mtx[mtx_ndx(miga[1],j,k)];
431 <                    float               mc = miga[2]->mtx[mtx_ndx(miga[2],i,k)];
430 >                    float               mb = mtx_coef(miga[1],j,k);
431 >                    float               mc = mtx_coef(miga[2],i,k);
432                      const RBFVAL        *rbf2k;
433                      double              rad2k;
464 –                    FVECT               vout;
434                      int                 pos[2];
435 <                    if ((mb <= FTINY) | (mc <= FTINY))
435 >                    if ((mb <= cthresh) & (mc <= cthresh))
436                          continue;
437                      rbf2k = &miga[2]->rbfv[1]->rbfa[k];
438                      rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact);
439                      rad2k = R2ANG(rbf2k->crad);
440                      rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
441                      ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
442 <                    geodesic(vout, v1, v2, t, GEOD_REL);
443 <                    pos_from_vec(pos, vout);
442 >                    geodesic(v2, v1, v2, t, GEOD_REL);
443 >                    pos_from_vec(pos, v2);
444                      rbf->rbfa[n].gx = pos[0];
445                      rbf->rbfa[n].gy = pos[1];
446                      ++n;

Diff Legend

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