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.2 by greg, Sat Oct 20 07:02:00 2012 UTC vs.
Revision 2.15 by greg, Wed Oct 23 03:41:39 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 +        if (cos_a >= 1.)                /* handles rounding error */
98 +                return(1);
99 +
100 +        cos_b = DOT(ej->rbfv[1]->invec, ivec);
101 +        if (cos_b <= 0)
102 +                return(0);
103 +        if (cos_b >= 1.)
104 +                return(1);
105 +
106 +        cos_aplusb = cos_a*cos_b - sqrt((1.-cos_a*cos_a)*(1.-cos_b*cos_b));
107 +        if (cos_aplusb <= 0)
108 +                return(0);
109 +
110 +        cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec);
111 +
112 +        return(cos_c - cos_aplusb < .001);
113 + }
114 +
115 + /* Determine if we are inside the given triangle */
116 + static int
117 + in_tri(const RBFNODE *v1, const RBFNODE *v2, const RBFNODE *v3, const FVECT p)
118 + {
119 +        FVECT   vc;
120 +        int     sgn1, sgn2, sgn3;
121 +                                        /* signed volume test */
122 +        VCROSS(vc, v1->invec, v2->invec);
123 +        sgn1 = (DOT(p, vc) > 0);
124 +        VCROSS(vc, v2->invec, v3->invec);
125 +        sgn2 = (DOT(p, vc) > 0);
126 +        if (sgn1 != sgn2)
127 +                return(0);
128 +        VCROSS(vc, v3->invec, v1->invec);
129 +        sgn3 = (DOT(p, vc) > 0);
130 +        return(sgn2 == sgn3);
131 + }
132 +
133 + /* Test (and set) bitmap for edge */
134 + static int
135 + check_edge(unsigned char *emap, int nedges, const MIGRATION *mig, int mark)
136 + {
137 +        int     ejndx, bit2check;
138 +
139 +        if (mig->rbfv[0]->ord > mig->rbfv[1]->ord)
140 +                ejndx = mig->rbfv[1]->ord + (nedges-1)*mig->rbfv[0]->ord;
141 +        else
142 +                ejndx = mig->rbfv[0]->ord + (nedges-1)*mig->rbfv[1]->ord;
143 +
144 +        bit2check = 1<<(ejndx&07);
145 +
146 +        if (emap[ejndx>>3] & bit2check)
147 +                return(0);
148 +        if (mark)
149 +                emap[ejndx>>3] |= bit2check;
150 +        return(1);
151 + }
152 +
153 + /* Compute intersection with the given position over remaining mesh */
154 + static int
155 + in_mesh(MIGRATION *miga[3], unsigned char *emap, int nedges,
156 +                        const FVECT ivec, MIGRATION *mig)
157 + {
158 +        RBFNODE         *tv[2];
159 +        MIGRATION       *sej[2], *dej[2];
160 +        int             i;
161 +                                                /* check visitation record */
162 +        if (!check_edge(emap, nedges, mig, 1))
163 +                return(0);
164 +        if (on_edge(mig, ivec)) {
165 +                miga[0] = mig;                  /* close enough to edge */
166 +                return(1);
167 +        }
168 +        if (!get_triangles(tv, mig))            /* do triangles either side? */
169 +                return(0);
170 +        for (i = 2; i--; ) {                    /* identify edges to check */
171 +                MIGRATION       *ej;
172 +                sej[i] = dej[i] = NULL;
173 +                if (tv[i] == NULL)
174 +                        continue;
175 +                for (ej = tv[i]->ejl; ej != NULL; ej = nextedge(tv[i],ej)) {
176 +                        RBFNODE *rbfop = opp_rbf(tv[i],ej);
177 +                        if (rbfop == mig->rbfv[0]) {
178 +                                if (check_edge(emap, nedges, ej, 0))
179 +                                        sej[i] = ej;
180 +                        } else if (rbfop == mig->rbfv[1]) {
181 +                                if (check_edge(emap, nedges, ej, 0))
182 +                                        dej[i] = ej;
183 +                        }
184 +                }
185 +        }
186 +        for (i = 2; i--; ) {                    /* check triangles just once */
187 +                if (sej[i] != NULL && in_mesh(miga, emap, nedges, ivec, sej[i]))
188 +                        return(1);
189 +                if (dej[i] != NULL && in_mesh(miga, emap, nedges, ivec, dej[i]))
190 +                        return(1);
191 +                if ((sej[i] == NULL) | (dej[i] == NULL))
192 +                        continue;
193 +                if (in_tri(mig->rbfv[0], mig->rbfv[1], tv[i], ivec)) {
194 +                        miga[0] = mig;
195 +                        miga[1] = sej[i];
196 +                        miga[2] = dej[i];
197 +                        return(1);
198 +                }
199 +        }
200 +        return(0);                              /* not near this edge */
201 + }
202 +
203   /* Find edge(s) for interpolating the given vector, applying symmetry */
204   int
205   get_interp(MIGRATION *miga[3], FVECT invec)
206   {
207          miga[0] = miga[1] = miga[2] = NULL;
208          if (single_plane_incident) {            /* isotropic BSDF? */
209 <                RBFNODE *rbf;                   /* find edge we're on */
210 <                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
211 <                        if (input_orient*rbf->invec[2] < input_orient*invec[2])
212 <                                break;
213 <                        if (rbf->next != NULL &&
233 <                                        input_orient*rbf->next->invec[2] <
209 >            RBFNODE     *rbf;                   /* find edge we're on */
210 >            for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
211 >                if (input_orient*rbf->invec[2] < input_orient*invec[2])
212 >                        break;
213 >                if (rbf->next != NULL && input_orient*rbf->next->invec[2] <
214                                                          input_orient*invec[2]) {
215 <                                for (miga[0] = rbf->ejl; miga[0] != NULL;
216 <                                                miga[0] = nextedge(rbf,miga[0]))
217 <                                        if (opp_rbf(rbf,miga[0]) == rbf->next)
218 <                                                return(0);
219 <                                break;
215 >                    for (miga[0] = rbf->ejl; miga[0] != NULL;
216 >                                        miga[0] = nextedge(rbf,miga[0]))
217 >                        if (opp_rbf(rbf,miga[0]) == rbf->next) {
218 >                                double  nf = 1. - rbf->invec[2]*rbf->invec[2];
219 >                                if (nf > FTINY) {       /* rotate to match */
220 >                                        nf = sqrt((1.-invec[2]*invec[2])/nf);
221 >                                        invec[0] = nf*rbf->invec[0];
222 >                                        invec[1] = nf*rbf->invec[1];
223 >                                }
224 >                                return(0);
225                          }
226 +                    break;
227                  }
228 <                return(-1);                     /* outside range! */
228 >            }
229 >            return(-1);                         /* outside range! */
230          }
231          {                                       /* else use triangle mesh */
232 <                const int       sym = use_symmetry(invec);
233 <                unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
234 <                int             pstart[2];
235 <                RBFNODE         *vother;
236 <                MIGRATION       *ej;
237 <                int             i;
238 <
239 <                pos_from_vec(pstart, invec);
240 <                memset(floodmap, 0, sizeof(floodmap));
241 <                                                /* call flooding function */
242 <                if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
243 <                        return(-1);             /* outside mesh */
244 <                if ((miga[0] == NULL) | (miga[2] == NULL))
245 <                        return(-1);             /* should never happen */
246 <                if (miga[1] == NULL)
260 <                        return(sym);            /* on edge */
261 <                                                /* verify triangle */
262 <                if (!order_triangle(miga)) {
232 >                int             sym = use_symmetry(invec);
233 >                int             nedges = 0;
234 >                MIGRATION       *mep;
235 >                unsigned char   *emap;
236 >                                                /* clear visitation map */
237 >                for (mep = mig_list; mep != NULL; mep = mep->next)
238 >                        ++nedges;
239 >                emap = (unsigned char *)calloc((nedges*(nedges-1) + 7)>>3, 1);
240 >                if (emap == NULL) {
241 >                        fprintf(stderr, "%s: Out of memory in get_interp()\n",
242 >                                        progname);
243 >                        exit(1);
244 >                }
245 >                                                /* identify intersection  */
246 >                if (!in_mesh(miga, emap, nedges, invec, mig_list)) {
247   #ifdef DEBUG
248 <                        fputs("Munged triangle in get_interp()\n", stderr);
248 >                        fprintf(stderr,
249 >                        "Incident angle (%.1f,%.1f) deg. outside mesh\n",
250 >                                        get_theta180(invec), get_phi360(invec));
251   #endif
252 <                        vother = NULL;          /* find triangle from edge */
253 <                        for (i = 3; i--; ) {
254 <                            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 */
252 >                        sym = -1;               /* outside mesh */
253 >                } else if (miga[1] != NULL &&
254 >                                (miga[2] == NULL || !order_triangle(miga))) {
255   #ifdef DEBUG
256 <                                fputs("No triangle in get_interp()\n", stderr);
256 >                        fputs("Munged triangle in get_interp()\n", stderr);
257   #endif
258 <                                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 <                        }
258 >                        sym = -1;
259                  }
260 +                free(emap);
261                  return(sym);                    /* return in standard order */
262          }
263   }
264  
265   /* Advect and allocate new RBF along edge */
266   static RBFNODE *
267 < e_advect_rbf(const MIGRATION *mig, const FVECT invec)
267 > e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
268   {
269 +        double          cthresh = FTINY;
270          RBFNODE         *rbf;
271          int             n, i, j;
272          double          t, full_dist;
273                                                  /* get relative position */
274 <        t = acos(DOT(invec, mig->rbfv[0]->invec));
275 <        if (t < M_PI/GRIDRES) {                 /* near first DSF */
274 >        t = Acos(DOT(invec, mig->rbfv[0]->invec));
275 >        if (t < M_PI/grid_res) {                /* near first DSF */
276                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
277                  rbf = (RBFNODE *)malloc(n);
278                  if (rbf == NULL)
279                          goto memerr;
280                  memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
281 +                rbf->next = NULL; rbf->ejl = NULL;
282                  return(rbf);
283          }
284          full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
285 <        if (t > full_dist-M_PI/GRIDRES) {       /* near second DSF */
285 >        if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
286                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
287                  rbf = (RBFNODE *)malloc(n);
288                  if (rbf == NULL)
289                          goto memerr;
290                  memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
291 +                rbf->next = NULL; rbf->ejl = NULL;
292                  return(rbf);
293          }
294 <        t /= full_dist;
294 >        t /= full_dist;
295 > tryagain:
296          n = 0;                                  /* count migrating particles */
297          for (i = 0; i < mtx_nrows(mig); i++)
298              for (j = 0; j < mtx_ncols(mig); j++)
299 <                n += (mtx_coef(mig,i,j) > FTINY);
299 >                n += (mtx_coef(mig,i,j) > cthresh);
300 >                                                /* are we over our limit? */
301 >        if ((lobe_lim > 0) & (n > lobe_lim)) {
302 >                cthresh = cthresh*2. + 10.*FTINY;
303 >                goto tryagain;
304 >        }
305   #ifdef DEBUG
306          fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
307                          mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
# Line 350 | Line 322 | e_advect_rbf(const MIGRATION *mig, const FVECT invec)
322              float               mv;
323              ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
324              for (j = 0; j < mtx_ncols(mig); j++)
325 <                if ((mv = mtx_coef(mig,i,j)) > FTINY) {
325 >                if ((mv = mtx_coef(mig,i,j)) > cthresh) {
326                          const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
327 <                        double          rad1 = R2ANG(rbf1j->crad);
327 >                        double          rad2;
328                          FVECT           v;
329                          int             pos[2];
330 <                        rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal;
331 <                        rbf->rbfa[n].crad = ANG2R(sqrt(rad0*rad0*(1.-t) +
332 <                                                        rad1*rad1*t));
330 >                        rad2 = R2ANG(rbf1j->crad);
331 >                        rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
332 >                        rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
333 >                                                rad0*rad0/rad2;
334 >                        rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
335                          ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
336                          geodesic(v, v0, v, t, GEOD_REL);
337                          pos_from_vec(pos, v);
# Line 376 | Line 350 | memerr:
350  
351   /* Partially advect between recorded incident angles and allocate new RBF */
352   RBFNODE *
353 < advect_rbf(const FVECT invec)
353 > advect_rbf(const FVECT invec, int lobe_lim)
354   {
355 +        double          cthresh = FTINY;
356          FVECT           sivec;
357          MIGRATION       *miga[3];
358          RBFNODE         *rbf;
# Line 392 | Line 367 | advect_rbf(const FVECT invec)
367          if (sym < 0)                            /* can't interpolate? */
368                  return(NULL);
369          if (miga[1] == NULL) {                  /* advect along edge? */
370 <                rbf = e_advect_rbf(miga[0], sivec);
371 <                rev_rbf_symmetry(rbf, sym);
370 >                rbf = e_advect_rbf(miga[0], sivec, lobe_lim);
371 >                if (single_plane_incident)
372 >                        rotate_rbf(rbf, invec);
373 >                else
374 >                        rev_rbf_symmetry(rbf, sym);
375                  return(rbf);
376          }
377   #ifdef DEBUG
# Line 415 | Line 393 | advect_rbf(const FVECT invec)
393          geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
394                          s, GEOD_REL);
395          t = acos(DOT(v1,sivec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
396 + tryagain:
397          n = 0;                                  /* count migrating particles */
398          for (i = 0; i < mtx_nrows(miga[0]); i++)
399              for (j = 0; j < mtx_ncols(miga[0]); j++)
400 <                for (k = (mtx_coef(miga[0],i,j) > FTINY) *
400 >                for (k = (mtx_coef(miga[0],i,j) > cthresh) *
401                                          mtx_ncols(miga[2]); k--; )
402 <                        n += (mtx_coef(miga[2],i,k) > FTINY &&
403 <                                mtx_coef(miga[1],j,k) > FTINY);
402 >                        n += (mtx_coef(miga[2],i,k) > cthresh ||
403 >                                mtx_coef(miga[1],j,k) > cthresh);
404 >                                                /* are we over our limit? */
405 >        if ((lobe_lim > 0) & (n > lobe_lim)) {
406 >                cthresh = cthresh*2. + 10.*FTINY;
407 >                goto tryagain;
408 >        }
409   #ifdef DEBUG
410          fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
411                          miga[0]->rbfv[0]->nrbf, miga[0]->rbfv[1]->nrbf,
# Line 448 | Line 432 | advect_rbf(const FVECT invec)
432              for (j = 0; j < mtx_ncols(miga[0]); j++) {
433                  const float     ma = mtx_coef(miga[0],i,j);
434                  const RBFVAL    *rbf1j;
435 <                double          rad1j, srad2;
436 <                if (ma <= FTINY)
435 >                double          srad2;
436 >                if (ma <= cthresh)
437                          continue;
438                  rbf1j = &miga[0]->rbfv[1]->rbfa[j];
439 <                rad1j = R2ANG(rbf1j->crad);
440 <                srad2 = (1.-s)*(1.-t)*rad0i*rad0i + s*(1.-t)*rad1j*rad1j;
439 >                srad2 = R2ANG(rbf1j->crad);
440 >                srad2 = (1.-s)*(1.-t)*rad0i*rad0i + s*(1.-t)*srad2*srad2;
441                  ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
442                  geodesic(v1, v0, v1, s, GEOD_REL);
443                  for (k = 0; k < mtx_ncols(miga[2]); k++) {
444                      float               mb = mtx_coef(miga[1],j,k);
445                      float               mc = mtx_coef(miga[2],i,k);
446                      const RBFVAL        *rbf2k;
447 <                    double              rad2k;
464 <                    FVECT               vout;
447 >                    double              rad2;
448                      int                 pos[2];
449 <                    if ((mb <= FTINY) | (mc <= FTINY))
449 >                    if ((mb <= cthresh) & (mc <= cthresh))
450                          continue;
451                      rbf2k = &miga[2]->rbfv[1]->rbfa[k];
452 <                    rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact);
453 <                    rad2k = R2ANG(rbf2k->crad);
454 <                    rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
452 >                    rad2 = R2ANG(rbf2k->crad);
453 >                    rad2 = srad2 + t*rad2*rad2;
454 >                    rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact) *
455 >                                        rad0i*rad0i/rad2;
456 >                    rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
457                      ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
458 <                    geodesic(vout, v1, v2, t, GEOD_REL);
459 <                    pos_from_vec(pos, vout);
458 >                    geodesic(v2, v1, v2, t, GEOD_REL);
459 >                    pos_from_vec(pos, v2);
460                      rbf->rbfa[n].gx = pos[0];
461                      rbf->rbfa[n].gy = pos[1];
462                      ++n;

Diff Legend

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