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

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

Diff Legend

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