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

Comparing ray/src/cv/bsdfrbf.c (file contents):
Revision 2.11 by greg, Sat Oct 19 00:11:50 2013 UTC vs.
Revision 2.12 by greg, Mon Oct 21 18:33:15 2013 UTC

# Line 14 | Line 14 | static const char RCSid[] = "$Id$";
14   #include <math.h>
15   #include "bsdfrep.h"
16  
17 < #ifndef MINRSCA
18 < #define MINRSCA         1.0             /* minimum radius scaling factor */
17 > #ifndef RSCA
18 > #define RSCA            2.2             /* radius scaling factor (empirical) */
19   #endif
20 < #ifndef MAXRSCA
21 < #define MAXRSCA         2.7             /* maximum radius scaling factor */
20 > #ifndef SMOOTH_MSE
21 > #define SMOOTH_MSE      5e-5            /* acceptable mean squared error */
22   #endif
23 < #ifndef VARTHRESH
24 < #define VARTHRESH       0.0015          /* culling variance threshold */
23 > #ifndef SMOOTH_MSER
24 > #define SMOOTH_MSER     0.07            /* acceptable relative MSE */
25   #endif
26 < #ifndef DIFFMAX2
27 < #define DIFFMAX2        (16.*VARTHRESH) /* maximum ignored sample variance */
28 < #endif
29 < #ifndef MAXFRAC
30 < #define MAXFRAC         0.5             /* maximum contribution to neighbor */
31 < #endif
32 < #ifndef NNEIGH
33 < #define NNEIGH          10              /* number of neighbors to consider */
34 < #endif
26 > #define MAX_RAD         (GRIDRES/8)     /* maximum lobe radius */
27 >
28 > #define RBFALLOCB       10              /* RBF allocation block size */
29 >
30                                  /* our loaded grid for this incident angle */
31   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
32  
# Line 78 | Line 73 | add_bsdf_data(double theta_out, double phi_out, double
73          dsf_grid[pos[0]][pos[1]].nval++;
74   }
75  
81 – /* Compute radii for non-empty bins */
82 – /* (distance to furthest empty bin for which non-empty test bin is closest) */
83 – static void
84 – compute_radii(void)
85 – {
86 –        const int       cradmin = ANG2R(.5*M_PI/GRIDRES);
87 –        unsigned int    fill_grid[GRIDRES][GRIDRES];
88 –        unsigned short  fill_cnt[GRIDRES][GRIDRES];
89 –        FVECT           ovec0, ovec1;
90 –        double          ang2, lastang2;
91 –        int             r, i, j, jn, ii, jj, inear, jnear;
92 –
93 –        r = GRIDRES/2;                          /* proceed in zig-zag */
94 –        for (i = 0; i < GRIDRES; i++)
95 –            for (jn = 0; jn < GRIDRES; jn++) {
96 –                j = (i&1) ? jn : GRIDRES-1-jn;
97 –                if (dsf_grid[i][j].nval)        /* find empty grid pos. */
98 –                        continue;
99 –                ovec_from_pos(ovec0, i, j);
100 –                inear = jnear = -1;             /* find nearest non-empty */
101 –                lastang2 = M_PI*M_PI;
102 –                for (ii = i-r; ii <= i+r; ii++) {
103 –                    if (ii < 0) continue;
104 –                    if (ii >= GRIDRES) break;
105 –                    for (jj = j-r; jj <= j+r; jj++) {
106 –                        if (jj < 0) continue;
107 –                        if (jj >= GRIDRES) break;
108 –                        if (!dsf_grid[ii][jj].nval)
109 –                                continue;
110 –                        ovec_from_pos(ovec1, ii, jj);
111 –                        ang2 = 2. - 2.*DOT(ovec0,ovec1);
112 –                        if (ang2 >= lastang2)
113 –                                continue;
114 –                        lastang2 = ang2;
115 –                        inear = ii; jnear = jj;
116 –                    }
117 –                }
118 –                if (inear < 0) {
119 –                        fprintf(stderr,
120 –                                "%s: Could not find non-empty neighbor!\n",
121 –                                        progname);
122 –                        exit(1);
123 –                }
124 –                ang2 = sqrt(lastang2);
125 –                r = ANG2R(ang2);                /* record if > previous */
126 –                if (r > dsf_grid[inear][jnear].crad)
127 –                        dsf_grid[inear][jnear].crad = r;
128 –                                                /* next search radius */
129 –                r = ang2*(2.*GRIDRES/M_PI) + 3;
130 –            }
131 –        for (i = 0; i < GRIDRES; i++)           /* grow radii where uniform */
132 –            for (j = 0; j < GRIDRES; j++) {
133 –                double  midmean = 0.0;
134 –                int     nsum = 0;
135 –                if (!dsf_grid[i][j].nval)
136 –                        continue;
137 –                r = 1;                          /* avg. immediate neighbors */
138 –                for (ii = i-r; ii <= i+r; ii++) {
139 –                    if (ii < 0) continue;
140 –                    if (ii >= GRIDRES) break;
141 –                    for (jj = j-r; jj <= j+r; jj++) {
142 –                        if (jj < 0) continue;
143 –                        if (jj >= GRIDRES) break;
144 –                        midmean += dsf_grid[ii][jj].vsum;
145 –                        nsum += dsf_grid[ii][jj].nval;
146 –                    }
147 –                }
148 –                midmean /= (double)nsum;
149 –                while (++r < GRIDRES) {         /* attempt to grow perimeter */
150 –                    double      diff2sum = 0.0;
151 –                    nsum = 0;
152 –                    for (ii = i-r; ii <= i+r; ii++) {
153 –                        int     jstep = 1;
154 –                        if (ii < 0) continue;
155 –                        if (ii >= GRIDRES) break;
156 –                        if ((i-r < ii) & (ii < i+r))
157 –                                jstep = r<<1;
158 –                        for (jj = j-r; jj <= j+r; jj += jstep) {
159 –                                double  d2;
160 –                                if (jj < 0) continue;
161 –                                if (jj >= GRIDRES) break;
162 –                                if (!dsf_grid[ii][jj].nval)
163 –                                        continue;
164 –                                d2 = midmean - dsf_grid[ii][jj].vsum /
165 –                                                (double)dsf_grid[ii][jj].nval;
166 –                                d2 *= d2;
167 –                                if (d2 > DIFFMAX2*midmean*midmean)
168 –                                        goto escape;
169 –                                diff2sum += d2;
170 –                                ++nsum;
171 –                        }
172 –                    }
173 –                    if (diff2sum > VARTHRESH*midmean*midmean*(double)nsum)
174 –                        break;
175 –                }
176 – escape:         --r;
177 –                r = ANG2R(r*(M_PI/MAXRSCA/GRIDRES));
178 –                if (r < cradmin)
179 –                        r = cradmin;
180 –                if (dsf_grid[i][j].crad < r)
181 –                        dsf_grid[i][j].crad = r;
182 –            }
183 –                                                /* blur radii over hemisphere */
184 –        memset(fill_grid, 0, sizeof(fill_grid));
185 –        memset(fill_cnt, 0, sizeof(fill_cnt));
186 –        for (i = 0; i < GRIDRES; i++)
187 –            for (j = 0; j < GRIDRES; j++) {
188 –                if (!dsf_grid[i][j].nval)
189 –                        continue;               /* not part of this */
190 –                r = R2ANG(dsf_grid[i][j].crad)*(2.*MAXRSCA*GRIDRES/M_PI);
191 –                for (ii = i-r; ii <= i+r; ii++) {
192 –                    if (ii < 0) continue;
193 –                    if (ii >= GRIDRES) break;
194 –                    for (jj = j-r; jj <= j+r; jj++) {
195 –                        if (jj < 0) continue;
196 –                        if (jj >= GRIDRES) break;
197 –                        if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
198 –                                continue;
199 –                        fill_grid[ii][jj] += dsf_grid[i][j].crad;
200 –                        fill_cnt[ii][jj]++;
201 –                    }
202 –                }
203 –            }
204 –                                                /* copy back blurred radii */
205 –        for (i = 0; i < GRIDRES; i++)
206 –            for (j = 0; j < GRIDRES; j++)
207 –                if (fill_cnt[i][j])
208 –                        dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
209 – }
210 –
211 – /* Radius comparison for qsort() */
212 – static int
213 – radius_cmp(const void *p1, const void *p2)
214 – {
215 –        return( (int)dsf_grid[0][*(const int *)p1].crad -
216 –                (int)dsf_grid[0][*(const int *)p2].crad );
217 – }
218 –
219 – /* Cull points for more uniform distribution, leave all nval 0 or 1 */
220 – static void
221 – cull_values(void)
222 – {
223 –        int     indx[GRIDRES*GRIDRES];
224 –        FVECT   ovec0, ovec1;
225 –        double  maxang, maxang2;
226 –        int     i, j, k, ii, jj, r;
227 –                                                /* sort by radius first */
228 –        for (k = GRIDRES*GRIDRES; k--; )
229 –                indx[k] = k;
230 –        qsort(indx, GRIDRES*GRIDRES, sizeof(int), &radius_cmp);
231 –                                                /* simple greedy algorithm */
232 –        for (k = GRIDRES*GRIDRES; k--; ) {
233 –                i = indx[k]/GRIDRES;            /* from biggest radius down */
234 –                j = indx[k] - i*GRIDRES;
235 –                if (!dsf_grid[i][j].nval)
236 –                        continue;
237 –                if (!dsf_grid[i][j].crad)
238 –                        break;          /* shouldn't happen */
239 –                ovec_from_pos(ovec0, i, j);
240 –                maxang = 2.*R2ANG(dsf_grid[i][j].crad);
241 –                                                /* clamp near horizon */
242 –                if (maxang > output_orient*ovec0[2])
243 –                        maxang = output_orient*ovec0[2];
244 –                r = maxang*(2.*GRIDRES/M_PI) + 1;
245 –                maxang2 = maxang*maxang;
246 –                for (ii = i-r; ii <= i+r; ii++) {
247 –                    if (ii < 0) continue;
248 –                    if (ii >= GRIDRES) break;
249 –                    for (jj = j-r; jj <= j+r; jj++) {
250 –                        if ((ii == i) & (jj == j))
251 –                                continue;       /* don't get self-absorbed */
252 –                        if (jj < 0) continue;
253 –                        if (jj >= GRIDRES) break;
254 –                        if (!dsf_grid[ii][jj].nval)
255 –                                continue;
256 –                        ovec_from_pos(ovec1, ii, jj);
257 –                        if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
258 –                                continue;
259 –                                                /* absorb sum */
260 –                        dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
261 –                        dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
262 –                                                /* keep value, though */
263 –                        dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
264 –                        dsf_grid[ii][jj].nval = 0;
265 –                    }
266 –                }
267 –        }
268 –                                                /* final averaging pass */
269 –        for (i = 0; i < GRIDRES; i++)
270 –            for (j = 0; j < GRIDRES; j++)
271 –                if (dsf_grid[i][j].nval > 1) {
272 –                        dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
273 –                        dsf_grid[i][j].nval = 1;
274 –                }
275 – }
276 –
76   /* Compute minimum BSDF from histogram (does not clear) */
77   static void
78   comp_bsdf_min()
# Line 295 | Line 94 | comp_bsdf_min()
94          bsdf_min = histval(i-1);
95   }
96  
97 < /* Find n nearest sub-sampled neighbors to the given grid position */
97 > /* Determine if the given region is empty of grid samples */
98   static int
99 < get_neighbors(int neigh[][2], int n, const int i, const int j)
99 > empty_region(int x0, int x1, int y0, int y1)
100   {
101 <        int     k = 0;
102 <        int     r;
103 <                                                /* search concentric squares */
104 <        for (r = 1; r < GRIDRES; r++) {
105 <                int     ii, jj;
106 <                for (ii = i-r; ii <= i+r; ii++) {
107 <                        int     jstep = 1;
108 <                        if (ii < 0) continue;
109 <                        if (ii >= GRIDRES) break;
110 <                        if ((i-r < ii) & (ii < i+r))
111 <                                jstep = r<<1;
112 <                        for (jj = j-r; jj <= j+r; jj += jstep) {
113 <                                if (jj < 0) continue;
114 <                                if (jj >= GRIDRES) break;
115 <                                if (dsf_grid[ii][jj].nval) {
116 <                                        neigh[k][0] = ii;
117 <                                        neigh[k][1] = jj;
118 <                                        if (++k >= n)
119 <                                                return(n);
120 <                                }
121 <                        }
101 >        int     x, y;
102 >
103 >        for (x = x0; x < x1; x++)
104 >            for (y = y0; y < y1; y++)
105 >                if (dsf_grid[x][y].nval)
106 >                        return(0);
107 >        return(1);
108 > }
109 >
110 > /* Determine if the given region is smooth enough to be a single lobe */
111 > static int
112 > smooth_region(int x0, int x1, int y0, int y1)
113 > {
114 >        RREAL   rMtx[3][3];
115 >        FVECT   xvec;
116 >        double  A, B, C, nvs, sqerr;
117 >        int     x, y, n;
118 >                                        /* compute planar regression */
119 >        memset(rMtx, 0, sizeof(rMtx));
120 >        memset(xvec, 0, sizeof(xvec));
121 >        for (x = x0; x < x1; x++)
122 >            for (y = y0; y < y1; y++)
123 >                if ((n = dsf_grid[x][y].nval) > 0) {
124 >                        double  z = dsf_grid[x][y].vsum;
125 >                        rMtx[0][0] += n*x*x;
126 >                        rMtx[0][1] += n*x*y;
127 >                        rMtx[0][2] += n*x;
128 >                        rMtx[1][1] += n*y*y;
129 >                        rMtx[1][2] += n*y;
130 >                        rMtx[2][2] += n;
131 >                        xvec[0] += x*z;
132 >                        xvec[1] += y*z;
133 >                        xvec[2] += z;
134                  }
135 +        rMtx[1][0] = rMtx[0][1];
136 +        rMtx[2][1] = rMtx[1][2];
137 +        nvs = rMtx[2][2];
138 +        if (SDinvXform(rMtx, rMtx) != SDEnone)
139 +                return(0);
140 +        A = DOT(rMtx[0], xvec);
141 +        B = DOT(rMtx[1], xvec);
142 +        C = DOT(rMtx[2], xvec);
143 +        sqerr = 0.0;                    /* compute mean squared error */
144 +        for (x = x0; x < x1; x++)
145 +            for (y = y0; y < y1; y++)
146 +                if ((n = dsf_grid[x][y].nval) > 0) {
147 +                        double  d = A*x + B*y + C - dsf_grid[x][y].vsum/n;
148 +                        sqerr += n*d*d;
149 +                }
150 +        if (sqerr <= nvs*SMOOTH_MSE)    /* below absolute MSE threshold? */
151 +                return(1);
152 +                                        /* below relative MSE threshold? */
153 +        return(sqerr*nvs <= xvec[2]*xvec[2]*SMOOTH_MSER);
154 + }
155 +
156 + /* Create new lobe based on integrated samples in region */
157 + static void
158 + create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
159 + {
160 +        double  vtot = 0.0;
161 +        int     nv = 0;
162 +        double  rad;
163 +        int     x, y;
164 +                                        /* compute average for region */
165 +        for (x = x0; x < x1; x++)
166 +            for (y = y0; y < y1; y++) {
167 +                vtot += dsf_grid[x][y].vsum;
168 +                nv += dsf_grid[x][y].nval;
169 +            }
170 +        if (!nv) {
171 +                fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
172 +                                progname);
173 +                exit(1);
174          }
175 <        return(k);
175 >                                        /* peak value based on integral */
176 >        vtot *= (x1-x0)*(y1-y0)*(2.*M_PI/GRIDRES/GRIDRES)/(double)nv;
177 >        rad = (RSCA/(double)GRIDRES)*(x1-x0);
178 >        rvp->peak =  vtot / ((2.*M_PI) * rad*rad);
179 >        rvp->crad = ANG2R(rad);
180 >        rvp->gx = (x0+x1)>>1;
181 >        rvp->gy = (y0+y1)>>1;
182   }
183  
184 < /* Adjust coded radius for the given grid position based on neighborhood */
184 > /* Recursive function to build radial basis function representation */
185   static int
186 < adj_coded_radius(const int i, const int j)
186 > build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, int y0, int y1)
187   {
188 <        const double    rad0 = R2ANG(dsf_grid[i][j].crad);
189 <        const double    minrad = MINRSCA * rad0;
190 <        double          currad = MAXRSCA * rad0;
191 <        int             neigh[NNEIGH][2];
192 <        int             n;
193 <        FVECT           our_dir;
194 <
195 <        ovec_from_pos(our_dir, i, j);
196 <        n = get_neighbors(neigh, NNEIGH, i, j);
197 <        while (n--) {
198 <                FVECT   their_dir;
199 <                double  max_ratio, rad_ok2;
200 <                                                /* check our value at neighbor */
201 <                ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
202 <                max_ratio = MAXFRAC * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
203 <                                / dsf_grid[i][j].vsum;
204 <                if (max_ratio >= 1)
205 <                        continue;
206 <                rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
207 <                if (rad_ok2 >= currad*currad)
208 <                        continue;               /* value fraction OK */
209 <                currad = sqrt(rad_ok2);         /* else reduce lobe radius */
210 <                if (currad <= minrad)           /* limit how small we'll go */
211 <                        return(ANG2R(minrad));
188 >        int     xmid = (x0+x1)>>1;
189 >        int     ymid = (y0+y1)>>1;
190 >        int     branched[4];
191 >        int     nadded, nleaves;
192 >                                        /* need to make this a leaf? */
193 >        if (empty_region(x0, xmid, y0, ymid) ||
194 >                        empty_region(xmid, x1, y0, ymid) ||
195 >                        empty_region(x0, xmid, ymid, y1) ||
196 >                        empty_region(xmid, x1, ymid, y1))
197 >                return(0);
198 >                                        /* add children (branches+leaves) */
199 >        if ((branched[0] = build_rbfrep(arp, np, x0, xmid, y0, ymid)) < 0)
200 >                return(-1);
201 >        if ((branched[1] = build_rbfrep(arp, np, xmid, x1, y0, ymid)) < 0)
202 >                return(-1);
203 >        if ((branched[2] = build_rbfrep(arp, np, x0, xmid, ymid, y1)) < 0)
204 >                return(-1);
205 >        if ((branched[3] = build_rbfrep(arp, np, xmid, x1, ymid, y1)) < 0)
206 >                return(-1);
207 >        nadded = branched[0] + branched[1] + branched[2] + branched[3];
208 >        nleaves = !branched[0] + !branched[1] + !branched[2] + !branched[3];
209 >        if (!nleaves)                   /* nothing but branches? */
210 >                return(nadded);
211 >                                        /* combine 4 leaves into 1? */
212 >        if (nleaves == 4 && x1-x0 < MAX_RAD && smooth_region(x0, x1, y0, y1))
213 >                return(0);
214 >                                        /* need more array space? */
215 >        if ((*np+nleaves-1)>>RBFALLOCB != (*np-1)>>RBFALLOCB) {
216 >                *arp = (RBFVAL *)realloc(*arp,
217 >                                sizeof(RBFVAL)*(*np+nleaves-1+(1<<RBFALLOCB)));
218 >                if (*arp == NULL)
219 >                        return(-1);
220          }
221 <        return(ANG2R(currad));                  /* encode selected radius */
221 >                                        /* create lobes for leaves */
222 >        if (!branched[0])
223 >                create_lobe(*arp+(*np)++, x0, xmid, y0, ymid);
224 >        if (!branched[1])
225 >                create_lobe(*arp+(*np)++, xmid, x1, y0, ymid);
226 >        if (!branched[2])
227 >                create_lobe(*arp+(*np)++, x0, xmid, ymid, y1);
228 >        if (!branched[3])
229 >                create_lobe(*arp+(*np)++, xmid, x1, ymid, y1);
230 >        nadded += nleaves;
231 >        return(nadded);
232   }
233  
234   /* Count up filled nodes and build RBF representation from current grid */
235   RBFNODE *
236 < make_rbfrep(void)
236 > make_rbfrep()
237   {
364 –        long    cradsum = 0, ocradsum = 0;
365 –        int     niter = 16;
366 –        double  lastVar, thisVar = 100.;
367 –        int     nn;
238          RBFNODE *newnode;
239 <        RBFVAL  *itera;
240 <        int     i, j;
371 <
372 < #ifdef DEBUG
373 < {
374 <        int     maxcnt = 0, nempty = 0;
375 <        long    cntsum = 0;
376 <        for (i = 0; i < GRIDRES; i++)
377 <            for (j = 0; j < GRIDRES; j++)
378 <                if (!dsf_grid[i][j].nval) {
379 <                        ++nempty;
380 <                } else {
381 <                        if (dsf_grid[i][j].nval > maxcnt)
382 <                                maxcnt = dsf_grid[i][j].nval;
383 <                        cntsum += dsf_grid[i][j].nval;
384 <                }
385 <        fprintf(stderr, "Average, maximum bin count: %d, %d (%.1f%% empty)\n",
386 <                        (int)(cntsum/((GRIDRES*GRIDRES)-nempty)), maxcnt,
387 <                        100./(GRIDRES*GRIDRES)*nempty);
388 < }
389 < #endif
390 <                                /* compute RBF radii */
391 <        compute_radii();
392 <                                /* coagulate lobes */
393 <        cull_values();
394 <        nn = 0;                 /* count selected bins */
395 <        for (i = 0; i < GRIDRES; i++)
396 <            for (j = 0; j < GRIDRES; j++)
397 <                nn += dsf_grid[i][j].nval;
239 >        RBFVAL  *rbfarr;
240 >        int     nn;
241                                  /* compute minimum BSDF */
242          comp_bsdf_min();
243 <                                /* allocate RBF array */
244 <        newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
243 >                                /* create RBF node list */
244 >        rbfarr = NULL; nn = 0;
245 >        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0)
246 >                goto memerr;
247 >                                /* (re)allocate RBF array */
248 >        newnode = (RBFNODE *)realloc(rbfarr,
249 >                        sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
250          if (newnode == NULL)
251                  goto memerr;
252 +                                /* copy computed lobes into RBF node */
253 +        memmove(newnode->rbfa, newnode, sizeof(RBFVAL)*nn);
254          newnode->ord = -1;
255          newnode->next = NULL;
256          newnode->ejl = NULL;
# Line 408 | Line 258 | make_rbfrep(void)
258          newnode->invec[0] = cos((M_PI/180.)*phi_in_deg)*newnode->invec[2];
259          newnode->invec[1] = sin((M_PI/180.)*phi_in_deg)*newnode->invec[2];
260          newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
261 <        newnode->vtotal = 0;
261 >        newnode->vtotal = .0;
262          newnode->nrbf = nn;
263 <        nn = 0;                 /* fill RBF array */
264 <        for (i = 0; i < GRIDRES; i++)
265 <            for (j = 0; j < GRIDRES; j++)
416 <                if (dsf_grid[i][j].nval) {
417 <                        newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
418 <                        ocradsum += dsf_grid[i][j].crad;
419 <                        cradsum +=
420 <                        newnode->rbfa[nn].crad = adj_coded_radius(i, j);
421 <                        newnode->rbfa[nn].gx = i;
422 <                        newnode->rbfa[nn].gy = j;
423 <                        ++nn;
424 <                }
263 >                                /* compute sum for normalization */
264 >        while (nn-- > 0)
265 >                newnode->vtotal += rbf_volume(&newnode->rbfa[nn]);
266   #ifdef DEBUG
267 <        fprintf(stderr,
427 <        "Average radius reduced from %.2f to %.2f degrees for %d lobes\n",
428 <                        180./M_PI*MAXRSCA*R2ANG(ocradsum/newnode->nrbf),
429 <                        180./M_PI*R2ANG(cradsum/newnode->nrbf), newnode->nrbf);
430 < #endif
431 <                                /* iterate to improve interpolation accuracy */
432 <        itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
433 <        if (itera == NULL)
434 <                goto memerr;
435 <        memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
436 <        do {
437 <                double  dsum = 0, dsum2 = 0;
438 <                nn = 0;
439 <                for (i = 0; i < GRIDRES; i++)
440 <                    for (j = 0; j < GRIDRES; j++)
441 <                        if (dsf_grid[i][j].nval) {
442 <                                FVECT   odir;
443 <                                double  corr;
444 <                                ovec_from_pos(odir, i, j);
445 <                                itera[nn++].peak *= corr =
446 <                                        dsf_grid[i][j].vsum /
447 <                                                eval_rbfrep(newnode, odir);
448 <                                dsum += 1. - corr;
449 <                                dsum2 += (1.-corr)*(1.-corr);
450 <                        }
451 <                memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
452 <                lastVar = thisVar;
453 <                thisVar = dsum2/(double)nn;
454 < #ifdef DEBUG
455 <                fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
456 <                                        100.*dsum/(double)nn,
457 <                                        100.*sqrt(thisVar));
458 < #endif
459 <        } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
460 <
461 <        free(itera);
462 <        nn = 0;                 /* compute sum for normalization */
463 <        while (nn < newnode->nrbf)
464 <                newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
465 < #ifdef DEBUG
267 >        fprintf(stderr, "Built RBF with %d lobes\n", newnode->nrbf);
268          fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
269                          get_theta180(newnode->invec), get_phi360(newnode->invec),
270                          newnode->vtotal);

Diff Legend

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