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.5 by greg, Fri Jun 28 23:18:51 2013 UTC vs.
Revision 2.28 by greg, Wed Feb 3 18:42:13 2016 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7   *      G. Ward
8   */
9  
10 + /****************************************************************
11 + 1) Collect samples into a grid using the Shirley-Chiu
12 +        angular mapping from a hemisphere to a square.
13 +
14 + 2) Compute an adaptive quadtree by subdividing the grid so that
15 +        each leaf node has at least one sample up to as many
16 +        samples as fit nicely on a plane to within a certain
17 +        MSE tolerance.
18 +
19 + 3) Place one Gaussian lobe at each leaf node in the quadtree,
20 +        sizing it to have a radius equal to the leaf size and
21 +        a volume equal to the energy in that node.
22 + *****************************************************************/
23 +
24   #define _USE_MATH_DEFINES
25   #include <stdio.h>
26   #include <stdlib.h>
# Line 15 | Line 29 | static const char RCSid[] = "$Id$";
29   #include "bsdfrep.h"
30  
31   #ifndef RSCA
32 < #define RSCA            2.7             /* radius scaling factor (empirical) */
32 > #define RSCA            2.0             /* radius scaling factor (empirical) */
33   #endif
34 <                                /* our loaded grid for this incident angle */
34 > #ifndef SMOOTH_MSE
35 > #define SMOOTH_MSE      5e-5            /* acceptable mean squared error */
36 > #endif
37 > #ifndef SMOOTH_MSER
38 > #define SMOOTH_MSER     0.03            /* acceptable relative MSE */
39 > #endif
40 > #define MAX_RAD         (GRIDRES/8)     /* maximum lobe radius */
41 >
42 > #define RBFALLOCB       10              /* RBF allocation block size */
43 >
44 >                                        /* loaded grid or comparison DSFs */
45   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
46 +                                        /* allocated chrominance sums if any */
47 + float                   (*spec_grid)[GRIDRES][GRIDRES];
48 + int                     nspec_grid = 0;
49  
50 + /* Set up visible spectrum sampling */
51 + void
52 + set_spectral_samples(int nspec)
53 + {
54 +        if (rbf_colorimetry == RBCunknown) {
55 +                if (nspec_grid > 0) {
56 +                        free(spec_grid);
57 +                        spec_grid = NULL;
58 +                        nspec_grid = 0;
59 +                }
60 +                if (nspec == 1) {
61 +                        rbf_colorimetry = RBCphotopic;
62 +                        return;
63 +                }
64 +                if (nspec == 3) {
65 +                        rbf_colorimetry = RBCtristimulus;
66 +                        spec_grid = (float (*)[GRIDRES][GRIDRES])calloc(
67 +                                        2*GRIDRES*GRIDRES, sizeof(float) );
68 +                        if (spec_grid == NULL)
69 +                                goto mem_error;
70 +                        nspec_grid = 2;
71 +                        return;
72 +                }
73 +                fprintf(stderr,
74 +                        "%s: only 1 or 3 spectral samples currently supported\n",
75 +                                progname);
76 +                exit(1);
77 +        }
78 +        if (nspec != nspec_grid+1) {
79 +                fprintf(stderr,
80 +                        "%s: number of spectral samples cannot be changed\n",
81 +                                progname);
82 +                exit(1);
83 +        }
84 +        return;
85 + mem_error:
86 +        fprintf(stderr, "%s: out of memory in set_spectral_samples()\n",
87 +                        progname);
88 +        exit(1);
89 + }
90 +
91   /* Start new DSF input grid */
92   void
93   new_bsdf_data(double new_theta, double new_phi)
# Line 27 | Line 95 | new_bsdf_data(double new_theta, double new_phi)
95          if (!new_input_direction(new_theta, new_phi))
96                  exit(1);
97          memset(dsf_grid, 0, sizeof(dsf_grid));
98 +        if (nspec_grid > 0)
99 +                memset(spec_grid, 0, sizeof(float)*GRIDRES*GRIDRES*nspec_grid);
100   }
101  
102   /* Add BSDF data point */
103   void
104 < add_bsdf_data(double theta_out, double phi_out, double val, int isDSF)
104 > add_bsdf_data(double theta_out, double phi_out, const double val[], int isDSF)
105   {
106          FVECT   ovec;
107 +        double  cfact, Yval;
108          int     pos[2];
109  
110 +        if (nspec_grid > 2) {
111 +                fprintf(stderr, "%s: unsupported color space\n", progname);
112 +                exit(1);
113 +        }
114          if (!output_orient)             /* check output orientation */
115                  output_orient = 1 - 2*(theta_out > 90.);
116          else if (output_orient > 0 ^ theta_out < 90.) {
117 <                fputs("Cannot handle output angles on both sides of surface\n",
118 <                                stderr);
117 >                fprintf(stderr,
118 >                "%s: cannot handle output angles on both sides of surface\n",
119 >                                progname);
120                  exit(1);
121          }
122          ovec[2] = sin((M_PI/180.)*theta_out);
123          ovec[0] = cos((M_PI/180.)*phi_out) * ovec[2];
124          ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
125          ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
126 +                                        /* BSDF to DSF correction */
127 +        cfact = isDSF ? 1. : ovec[2];
128  
129 <        if (!isDSF)
52 <                val *= ovec[2];         /* convert from BSDF to DSF */
53 <
129 >        Yval = cfact * val[rbf_colorimetry==RBCtristimulus];
130                                          /* update BSDF histogram */
131 <        if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
132 <                ++bsdf_hist[histndx(val/ovec[2])];
131 >        if (BSDF2SML*ovec[2] < Yval && Yval < BSDF2BIG*ovec[2])
132 >                ++bsdf_hist[histndx(Yval/ovec[2])];
133  
134          pos_from_vec(pos, ovec);
135  
136 <        dsf_grid[pos[0]][pos[1]].vsum += val;
137 <        dsf_grid[pos[0]][pos[1]].nval++;
136 >        dsf_grid[pos[0]][pos[1]].sum.v += Yval;
137 >        dsf_grid[pos[0]][pos[1]].sum.n++;
138 >                                        /* add in X and Z values */
139 >        if (rbf_colorimetry == RBCtristimulus) {
140 >                spec_grid[0][pos[0]][pos[1]] += cfact * val[0];
141 >                spec_grid[1][pos[0]][pos[1]] += cfact * val[2];
142 >        }
143   }
144  
145 < /* Compute radii for non-empty bins */
65 < /* (distance to furthest empty bin for which non-empty bin is the closest) */
145 > /* Compute minimum BSDF from histogram (does not clear) */
146   static void
67 – compute_radii(void)
68 – {
69 –        unsigned int    fill_grid[GRIDRES][GRIDRES];
70 –        unsigned short  fill_cnt[GRIDRES][GRIDRES];
71 –        FVECT           ovec0, ovec1;
72 –        double          ang2, lastang2;
73 –        int             r, i, j, jn, ii, jj, inear, jnear;
74 –
75 –        r = GRIDRES/2;                          /* proceed in zig-zag */
76 –        for (i = 0; i < GRIDRES; i++)
77 –            for (jn = 0; jn < GRIDRES; jn++) {
78 –                j = (i&1) ? jn : GRIDRES-1-jn;
79 –                if (dsf_grid[i][j].nval)        /* find empty grid pos. */
80 –                        continue;
81 –                ovec_from_pos(ovec0, i, j);
82 –                inear = jnear = -1;             /* find nearest non-empty */
83 –                lastang2 = M_PI*M_PI;
84 –                for (ii = i-r; ii <= i+r; ii++) {
85 –                    if (ii < 0) continue;
86 –                    if (ii >= GRIDRES) break;
87 –                    for (jj = j-r; jj <= j+r; jj++) {
88 –                        if (jj < 0) continue;
89 –                        if (jj >= GRIDRES) break;
90 –                        if (!dsf_grid[ii][jj].nval)
91 –                                continue;
92 –                        ovec_from_pos(ovec1, ii, jj);
93 –                        ang2 = 2. - 2.*DOT(ovec0,ovec1);
94 –                        if (ang2 >= lastang2)
95 –                                continue;
96 –                        lastang2 = ang2;
97 –                        inear = ii; jnear = jj;
98 –                    }
99 –                }
100 –                if (inear < 0) {
101 –                        fprintf(stderr,
102 –                                "%s: Could not find non-empty neighbor!\n",
103 –                                        progname);
104 –                        exit(1);
105 –                }
106 –                ang2 = sqrt(lastang2);
107 –                r = ANG2R(ang2);                /* record if > previous */
108 –                if (r > dsf_grid[inear][jnear].crad)
109 –                        dsf_grid[inear][jnear].crad = r;
110 –                                                /* next search radius */
111 –                r = ang2*(2.*GRIDRES/M_PI) + 3;
112 –            }
113 –                                                /* blur radii over hemisphere */
114 –        memset(fill_grid, 0, sizeof(fill_grid));
115 –        memset(fill_cnt, 0, sizeof(fill_cnt));
116 –        for (i = 0; i < GRIDRES; i++)
117 –            for (j = 0; j < GRIDRES; j++) {
118 –                if (!dsf_grid[i][j].crad)
119 –                        continue;               /* missing distance */
120 –                r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
121 –                for (ii = i-r; ii <= i+r; ii++) {
122 –                    if (ii < 0) continue;
123 –                    if (ii >= GRIDRES) break;
124 –                    for (jj = j-r; jj <= j+r; jj++) {
125 –                        if (jj < 0) continue;
126 –                        if (jj >= GRIDRES) break;
127 –                        if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
128 –                                continue;
129 –                        fill_grid[ii][jj] += dsf_grid[i][j].crad;
130 –                        fill_cnt[ii][jj]++;
131 –                    }
132 –                }
133 –            }
134 –                                                /* copy back blurred radii */
135 –        for (i = 0; i < GRIDRES; i++)
136 –            for (j = 0; j < GRIDRES; j++)
137 –                if (fill_cnt[i][j])
138 –                        dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
139 – }
140 –
141 – /* Cull points for more uniform distribution, leave all nval 0 or 1 */
142 – static void
143 – cull_values(void)
144 – {
145 –        FVECT   ovec0, ovec1;
146 –        double  maxang, maxang2;
147 –        int     i, j, ii, jj, r;
148 –                                                /* simple greedy algorithm */
149 –        for (i = 0; i < GRIDRES; i++)
150 –            for (j = 0; j < GRIDRES; j++) {
151 –                if (!dsf_grid[i][j].nval)
152 –                        continue;
153 –                if (!dsf_grid[i][j].crad)
154 –                        continue;               /* shouldn't happen */
155 –                ovec_from_pos(ovec0, i, j);
156 –                maxang = 2.*R2ANG(dsf_grid[i][j].crad);
157 –                if (maxang > ovec0[2])          /* clamp near horizon */
158 –                        maxang = ovec0[2];
159 –                r = maxang*(2.*GRIDRES/M_PI) + 1;
160 –                maxang2 = maxang*maxang;
161 –                for (ii = i-r; ii <= i+r; ii++) {
162 –                    if (ii < 0) continue;
163 –                    if (ii >= GRIDRES) break;
164 –                    for (jj = j-r; jj <= j+r; jj++) {
165 –                        if (jj < 0) continue;
166 –                        if (jj >= GRIDRES) break;
167 –                        if (!dsf_grid[ii][jj].nval)
168 –                                continue;
169 –                        if ((ii == i) & (jj == j))
170 –                                continue;       /* don't get self-absorbed */
171 –                        ovec_from_pos(ovec1, ii, jj);
172 –                        if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
173 –                                continue;
174 –                                                /* absorb sum */
175 –                        dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
176 –                        dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
177 –                                                /* keep value, though */
178 –                        dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
179 –                        dsf_grid[ii][jj].nval = 0;
180 –                    }
181 –                }
182 –            }
183 –                                                /* final averaging pass */
184 –        for (i = 0; i < GRIDRES; i++)
185 –            for (j = 0; j < GRIDRES; j++)
186 –                if (dsf_grid[i][j].nval > 1) {
187 –                        dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
188 –                        dsf_grid[i][j].nval = 1;
189 –                }
190 – }
191 –
192 – /* Compute minimum BSDF from histogram and clear it */
193 – static void
147   comp_bsdf_min()
148   {
149 <        int     cnt;
150 <        int     i, target;
149 >        unsigned long   cnt, target;
150 >        int             i;
151  
152          cnt = 0;
153          for (i = HISTLEN; i--; )
# Line 208 | Line 161 | comp_bsdf_min()
161          for (i = 0; cnt <= target; i++)
162                  cnt += bsdf_hist[i];
163          bsdf_min = histval(i-1);
211 –        memset(bsdf_hist, 0, sizeof(bsdf_hist));
164   }
165  
166 + /* Determine if the given region is empty of grid samples */
167 + static int
168 + empty_region(int x0, int x1, int y0, int y1)
169 + {
170 +        int     x, y;
171 +
172 +        for (x = x0; x < x1; x++)
173 +            for (y = y0; y < y1; y++)
174 +                if (dsf_grid[x][y].sum.n)
175 +                        return(0);
176 +        return(1);
177 + }
178 +
179 + /* Determine if the given region is smooth enough to be a single lobe */
180 + static int
181 + smooth_region(int x0, int x1, int y0, int y1)
182 + {
183 +        RREAL   rMtx[3][3];
184 +        FVECT   xvec;
185 +        double  A, B, C, nvs, sqerr;
186 +        int     x, y, n;
187 +                                        /* compute planar regression */
188 +        memset(rMtx, 0, sizeof(rMtx));
189 +        memset(xvec, 0, sizeof(xvec));
190 +        for (x = x0; x < x1; x++)
191 +            for (y = y0; y < y1; y++)
192 +                if ((n = dsf_grid[x][y].sum.n) > 0) {
193 +                        double  z = dsf_grid[x][y].sum.v;
194 +                        rMtx[0][0] += x*x*(double)n;
195 +                        rMtx[0][1] += x*y*(double)n;
196 +                        rMtx[0][2] += x*(double)n;
197 +                        rMtx[1][1] += y*y*(double)n;
198 +                        rMtx[1][2] += y*(double)n;
199 +                        rMtx[2][2] += (double)n;
200 +                        xvec[0] += x*z;
201 +                        xvec[1] += y*z;
202 +                        xvec[2] += z;
203 +                }
204 +        rMtx[1][0] = rMtx[0][1];
205 +        rMtx[2][0] = rMtx[0][2];
206 +        rMtx[2][1] = rMtx[1][2];
207 +        nvs = rMtx[2][2];
208 +        if (SDinvXform(rMtx, rMtx) != SDEnone)
209 +                return(1);              /* colinear values */
210 +        A = DOT(rMtx[0], xvec);
211 +        B = DOT(rMtx[1], xvec);
212 +        C = DOT(rMtx[2], xvec);
213 +        sqerr = 0.0;                    /* compute mean squared error */
214 +        for (x = x0; x < x1; x++)
215 +            for (y = y0; y < y1; y++)
216 +                if ((n = dsf_grid[x][y].sum.n) > 0) {
217 +                        double  d = A*x + B*y + C - dsf_grid[x][y].sum.v/n;
218 +                        sqerr += n*d*d;
219 +                }
220 +        if (sqerr <= nvs*SMOOTH_MSE)    /* below absolute MSE threshold? */
221 +                return(1);
222 +                                        /* OR below relative MSE threshold? */
223 +        return(sqerr*nvs <= xvec[2]*xvec[2]*SMOOTH_MSER);
224 + }
225 +
226 + /* Create new lobe based on integrated samples in region */
227 + static int
228 + create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
229 + {
230 +        double  vtot = 0.0;
231 +        double  CIEXtot = 0.0, CIEZtot = 0.0;
232 +        int     nv = 0;
233 +        double  wxsum = 0.0, wysum = 0.0, wtsum = 0.0;
234 +        double  rad;
235 +        int     x, y;
236 +                                        /* compute average for region */
237 +        for (x = x0; x < x1; x++)
238 +            for (y = y0; y < y1; y++)
239 +                if (dsf_grid[x][y].sum.n) {
240 +                    const double        v = dsf_grid[x][y].sum.v;
241 +                    const int           n = dsf_grid[x][y].sum.n;
242 +
243 +                    if (v > 0) {
244 +                        const double    wt = v / (double)n;
245 +                        wxsum += wt * x;
246 +                        wysum += wt * y;
247 +                        wtsum += wt;
248 +                    }
249 +                    vtot += v;
250 +                    nv += n;
251 +                    if (rbf_colorimetry == RBCtristimulus) {
252 +                        CIEXtot += spec_grid[0][x][y];
253 +                        CIEZtot += spec_grid[1][x][y];
254 +                    }
255 +                }
256 +        if (!nv) {
257 +                fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
258 +                                progname);
259 +                exit(1);
260 +        }
261 +        if (vtot <= 0)                  /* only create positive lobes */
262 +                return(0);
263 +                                        /* assign color */
264 +        if (rbf_colorimetry == RBCtristimulus) {
265 +                const double    df = 1.0 / (CIEXtot + vtot + CIEZtot);
266 +                C_COLOR         cclr;
267 +                c_cset(&cclr, CIEXtot*df, vtot*df);
268 +                rvp->chroma = c_encodeChroma(&cclr);
269 +        } else
270 +                rvp->chroma = c_dfchroma;
271 +                                        /* peak value based on integral */
272 +        vtot *= (x1-x0)*(y1-y0)*(2.*M_PI/GRIDRES/GRIDRES)/(double)nv;
273 +        rad = (RSCA/(double)GRIDRES)*(x1-x0);
274 +        rvp->peak =  vtot / ((2.*M_PI) * rad*rad);
275 +        rvp->crad = ANG2R(rad);         /* put peak at centroid */
276 +        rvp->gx = (int)(wxsum/wtsum + .5);
277 +        rvp->gy = (int)(wysum/wtsum + .5);
278 +        return(1);
279 + }
280 +
281 + /* Recursive function to build radial basis function representation */
282 + static int
283 + build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, int y0, int y1)
284 + {
285 +        int     xmid = (x0+x1)>>1;
286 +        int     ymid = (y0+y1)>>1;
287 +        int     branched[4];
288 +        int     nadded, nleaves;
289 +                                        /* need to make this a leaf? */
290 +        if (empty_region(x0, xmid, y0, ymid) ||
291 +                        empty_region(xmid, x1, y0, ymid) ||
292 +                        empty_region(x0, xmid, ymid, y1) ||
293 +                        empty_region(xmid, x1, ymid, y1))
294 +                return(0);
295 +                                        /* add children (branches+leaves) */
296 +        if ((branched[0] = build_rbfrep(arp, np, x0, xmid, y0, ymid)) < 0)
297 +                return(-1);
298 +        if ((branched[1] = build_rbfrep(arp, np, xmid, x1, y0, ymid)) < 0)
299 +                return(-1);
300 +        if ((branched[2] = build_rbfrep(arp, np, x0, xmid, ymid, y1)) < 0)
301 +                return(-1);
302 +        if ((branched[3] = build_rbfrep(arp, np, xmid, x1, ymid, y1)) < 0)
303 +                return(-1);
304 +        nadded = branched[0] + branched[1] + branched[2] + branched[3];
305 +        nleaves = !branched[0] + !branched[1] + !branched[2] + !branched[3];
306 +        if (!nleaves)                   /* nothing but branches? */
307 +                return(nadded);
308 +                                        /* combine 4 leaves into 1? */
309 +        if ((nleaves == 4) & (x1-x0 <= MAX_RAD) &&
310 +                        smooth_region(x0, x1, y0, y1))
311 +                return(0);
312 +                                        /* need more array space? */
313 +        if ((*np+nleaves-1)>>RBFALLOCB != (*np-1)>>RBFALLOCB) {
314 +                *arp = (RBFVAL *)realloc(*arp,
315 +                                sizeof(RBFVAL)*(*np+nleaves-1+(1<<RBFALLOCB)));
316 +                if (*arp == NULL)
317 +                        return(-1);
318 +        }
319 +                                        /* create lobes for leaves */
320 +        if (!branched[0] && create_lobe(*arp+*np, x0, xmid, y0, ymid)) {
321 +                ++(*np); ++nadded;
322 +        }
323 +        if (!branched[1] && create_lobe(*arp+*np, xmid, x1, y0, ymid)) {
324 +                ++(*np); ++nadded;
325 +        }
326 +        if (!branched[2] && create_lobe(*arp+*np, x0, xmid, ymid, y1)) {
327 +                ++(*np); ++nadded;
328 +        }
329 +        if (!branched[3] && create_lobe(*arp+*np, xmid, x1, ymid, y1)) {
330 +                ++(*np); ++nadded;
331 +        }
332 +        return(nadded);
333 + }
334 +
335   /* Count up filled nodes and build RBF representation from current grid */
336   RBFNODE *
337 < make_rbfrep(void)
337 > make_rbfrep()
338   {
218 –        int     niter = 16;
219 –        double  lastVar, thisVar = 100.;
220 –        int     nn;
339          RBFNODE *newnode;
340 <        RBFVAL  *itera;
341 <        int     i, j;
224 <                                /* compute RBF radii */
225 <        compute_radii();
226 <                                /* coagulate lobes */
227 <        cull_values();
228 <        nn = 0;                 /* count selected bins */
229 <        for (i = 0; i < GRIDRES; i++)
230 <            for (j = 0; j < GRIDRES; j++)
231 <                nn += dsf_grid[i][j].nval;
340 >        RBFVAL  *rbfarr;
341 >        int     nn;
342                                  /* compute minimum BSDF */
343          comp_bsdf_min();
344 <                                /* allocate RBF array */
345 <        newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
344 >                                /* create RBF node list */
345 >        rbfarr = NULL; nn = 0;
346 >        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0) {
347 >                if (nn)
348 >                        goto memerr;
349 >                fprintf(stderr,
350 >                        "%s: warning - skipping bad incidence (%.1f,%.1f)\n",
351 >                                progname, theta_in_deg, phi_in_deg);
352 >                return(NULL);
353 >        }
354 >                                /* (re)allocate RBF array */
355 >        newnode = (RBFNODE *)realloc(rbfarr,
356 >                        sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
357          if (newnode == NULL)
358                  goto memerr;
359 +                                /* copy computed lobes into RBF node */
360 +        memmove(newnode->rbfa, newnode, sizeof(RBFVAL)*nn);
361          newnode->ord = -1;
362          newnode->next = NULL;
363          newnode->ejl = NULL;
# Line 242 | Line 365 | make_rbfrep(void)
365          newnode->invec[0] = cos((M_PI/180.)*phi_in_deg)*newnode->invec[2];
366          newnode->invec[1] = sin((M_PI/180.)*phi_in_deg)*newnode->invec[2];
367          newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
368 <        newnode->vtotal = 0;
368 >        newnode->vtotal = .0;
369          newnode->nrbf = nn;
370 <        nn = 0;                 /* fill RBF array */
371 <        for (i = 0; i < GRIDRES; i++)
372 <            for (j = 0; j < GRIDRES; j++)
250 <                if (dsf_grid[i][j].nval) {
251 <                        newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
252 <                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
253 <                        newnode->rbfa[nn].gx = i;
254 <                        newnode->rbfa[nn].gy = j;
255 <                        ++nn;
256 <                }
257 <                                /* iterate to improve interpolation accuracy */
258 <        itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
259 <        if (itera == NULL)
260 <                goto memerr;
261 <        memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
262 <        do {
263 <                double  dsum = 0, dsum2 = 0;
264 <                nn = 0;
265 <                for (i = 0; i < GRIDRES; i++)
266 <                    for (j = 0; j < GRIDRES; j++)
267 <                        if (dsf_grid[i][j].nval) {
268 <                                FVECT   odir;
269 <                                double  corr;
270 <                                ovec_from_pos(odir, i, j);
271 <                                itera[nn++].peak *= corr =
272 <                                        dsf_grid[i][j].vsum /
273 <                                                eval_rbfrep(newnode, odir);
274 <                                dsum += 1. - corr;
275 <                                dsum2 += (1.-corr)*(1.-corr);
276 <                        }
277 <                memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
278 <                lastVar = thisVar;
279 <                thisVar = dsum2/(double)nn;
370 >                                /* compute sum for normalization */
371 >        while (nn-- > 0)
372 >                newnode->vtotal += rbf_volume(&newnode->rbfa[nn]);
373   #ifdef DEBUG
374 <                fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
282 <                                        100.*dsum/(double)nn,
283 <                                        100.*sqrt(thisVar));
284 < #endif
285 <        } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
286 <
287 <        free(itera);
288 <        nn = 0;                 /* compute sum for normalization */
289 <        while (nn < newnode->nrbf)
290 <                newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
291 < #ifdef DEBUG
374 >        fprintf(stderr, "Built RBF with %d lobes\n", newnode->nrbf);
375          fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
376                          get_theta180(newnode->invec), get_phi360(newnode->invec),
377                          newnode->vtotal);
378   #endif
379          insert_dsf(newnode);
297 –
380          return(newnode);
381   memerr:
382          fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);

Diff Legend

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