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.10 by greg, Fri Oct 18 02:49:30 2013 UTC vs.
Revision 2.26 by greg, Fri Aug 22 05:38:44 2014 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 14 | Line 28 | static const char RCSid[] = "$Id$";
28   #include <math.h>
29   #include "bsdfrep.h"
30  
31 < #ifndef MINRSCA
32 < #define MINRSCA         0.5             /* minimum radius scaling factor */
31 > #ifndef RSCA
32 > #define RSCA            2.0             /* radius scaling factor (empirical) */
33   #endif
34 < #ifndef MAXRSCA
35 < #define MAXRSCA         2.7             /* maximum radius scaling factor */
34 > #ifndef SMOOTH_MSE
35 > #define SMOOTH_MSE      5e-5            /* acceptable mean squared error */
36   #endif
37 < #ifndef DIFFTHRESH
38 < #define DIFFTHRESH      0.2             /* culling difference threshold */
37 > #ifndef SMOOTH_MSER
38 > #define SMOOTH_MSER     0.03            /* acceptable relative MSE */
39   #endif
40 < #ifndef MAXFRAC
41 < #define MAXFRAC         0.5             /* maximum contribution to neighbor */
42 < #endif
43 < #ifndef NNEIGH
44 < #define NNEIGH          10              /* number of neighbors to consider */
31 < #endif
32 <                                /* our loaded grid for this incident angle */
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  
47   /* Start new DSF input grid */
# Line 60 | Line 72 | add_bsdf_data(double theta_out, double phi_out, double
72          ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
73          ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
74  
75 <        if (val <= 0)                   /* truncate to zero */
76 <                val = 0;
65 <        else if (!isDSF)
66 <                val *= ovec[2];         /* convert from BSDF to DSF */
75 >        if (!isDSF)
76 >                val *= COSF(ovec[2]);   /* convert from BSDF to DSF */
77  
78                                          /* update BSDF histogram */
79          if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
# Line 71 | Line 81 | add_bsdf_data(double theta_out, double phi_out, double
81  
82          pos_from_vec(pos, ovec);
83  
84 <        dsf_grid[pos[0]][pos[1]].vsum += val;
85 <        dsf_grid[pos[0]][pos[1]].nval++;
84 >        dsf_grid[pos[0]][pos[1]].sum.v += val;
85 >        dsf_grid[pos[0]][pos[1]].sum.n++;
86   }
87  
88 < /* Check if the two DSF values are significantly different */
79 < static int
80 < big_diff(double ref, double tst)
81 < {
82 <        if (ref > 0) {
83 <                tst = tst/ref - 1.;
84 <                if (tst < 0) tst = -tst;
85 <        } else
86 <                tst *= 50.;
87 <        return(tst > DIFFTHRESH);
88 < }
89 <
90 < /* Compute radii for non-empty bins */
91 < /* (distance to furthest empty bin for which non-empty test bin is closest) */
88 > /* Compute minimum BSDF from histogram (does not clear) */
89   static void
93 – compute_radii(void)
94 – {
95 –        const int       cradmin = ANG2R(.5*M_PI/GRIDRES);
96 –        unsigned int    fill_grid[GRIDRES][GRIDRES];
97 –        unsigned short  fill_cnt[GRIDRES][GRIDRES];
98 –        FVECT           ovec0, ovec1;
99 –        double          ang2, lastang2;
100 –        int             r, i, j, jn, ii, jj, inear, jnear;
101 –
102 –        r = GRIDRES/2;                          /* proceed in zig-zag */
103 –        for (i = 0; i < GRIDRES; i++)
104 –            for (jn = 0; jn < GRIDRES; jn++) {
105 –                j = (i&1) ? jn : GRIDRES-1-jn;
106 –                if (dsf_grid[i][j].nval)        /* find empty grid pos. */
107 –                        continue;
108 –                ovec_from_pos(ovec0, i, j);
109 –                inear = jnear = -1;             /* find nearest non-empty */
110 –                lastang2 = M_PI*M_PI;
111 –                for (ii = i-r; ii <= i+r; ii++) {
112 –                    if (ii < 0) continue;
113 –                    if (ii >= GRIDRES) break;
114 –                    for (jj = j-r; jj <= j+r; jj++) {
115 –                        if (jj < 0) continue;
116 –                        if (jj >= GRIDRES) break;
117 –                        if (!dsf_grid[ii][jj].nval)
118 –                                continue;
119 –                        ovec_from_pos(ovec1, ii, jj);
120 –                        ang2 = 2. - 2.*DOT(ovec0,ovec1);
121 –                        if (ang2 >= lastang2)
122 –                                continue;
123 –                        lastang2 = ang2;
124 –                        inear = ii; jnear = jj;
125 –                    }
126 –                }
127 –                if (inear < 0) {
128 –                        fprintf(stderr,
129 –                                "%s: Could not find non-empty neighbor!\n",
130 –                                        progname);
131 –                        exit(1);
132 –                }
133 –                ang2 = sqrt(lastang2);
134 –                r = ANG2R(ang2);                /* record if > previous */
135 –                if (r > dsf_grid[inear][jnear].crad)
136 –                        dsf_grid[inear][jnear].crad = r;
137 –                                                /* next search radius */
138 –                r = ang2*(2.*GRIDRES/M_PI) + 3;
139 –            }
140 –        for (i = 0; i < GRIDRES; i++)           /* grow radii where uniform */
141 –            for (j = 0; j < GRIDRES; j++) {
142 –                double  midmean;
143 –                if (!dsf_grid[i][j].nval)
144 –                        continue;
145 –                midmean = dsf_grid[i][j].vsum / (double)dsf_grid[i][j].nval;
146 –                r = R2ANG(dsf_grid[i][j].crad)*(MAXRSCA*GRIDRES/M_PI);
147 –                while (++r < GRIDRES) {         /* attempt to grow perimeter */
148 –                    for (ii = i-r; ii <= i+r; ii++) {
149 –                        int     jstep = 1;
150 –                        if (ii < 0) continue;
151 –                        if (ii >= GRIDRES) break;
152 –                        if ((i-r < ii) & (ii < i+r))
153 –                                jstep = r<<1;
154 –                        for (jj = j-r; jj <= j+r; jj += jstep) {
155 –                                if (jj < 0) continue;
156 –                                if (jj >= GRIDRES) break;
157 –                                if (dsf_grid[ii][jj].nval && big_diff(midmean,
158 –                                                dsf_grid[ii][jj].vsum /
159 –                                                (double)dsf_grid[ii][jj].nval))
160 –                                        goto hit_diff;
161 –                        }
162 –                    }
163 –                }
164 – hit_diff:       --r;
165 –                dsf_grid[i][j].crad = ANG2R(r*(M_PI/MAXRSCA/GRIDRES));
166 –                if (dsf_grid[i][j].crad < cradmin)
167 –                        dsf_grid[i][j].crad = cradmin;
168 –            }
169 –                                                /* blur radii over hemisphere */
170 –        memset(fill_grid, 0, sizeof(fill_grid));
171 –        memset(fill_cnt, 0, sizeof(fill_cnt));
172 –        for (i = 0; i < GRIDRES; i++)
173 –            for (j = 0; j < GRIDRES; j++) {
174 –                if (!dsf_grid[i][j].nval)
175 –                        continue;               /* not part of this */
176 –                r = R2ANG(dsf_grid[i][j].crad)*(2.*MAXRSCA*GRIDRES/M_PI);
177 –                for (ii = i-r; ii <= i+r; ii++) {
178 –                    if (ii < 0) continue;
179 –                    if (ii >= GRIDRES) break;
180 –                    for (jj = j-r; jj <= j+r; jj++) {
181 –                        if (jj < 0) continue;
182 –                        if (jj >= GRIDRES) break;
183 –                        if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
184 –                                continue;
185 –                        fill_grid[ii][jj] += dsf_grid[i][j].crad;
186 –                        fill_cnt[ii][jj]++;
187 –                    }
188 –                }
189 –            }
190 –                                                /* copy back blurred radii */
191 –        for (i = 0; i < GRIDRES; i++)
192 –            for (j = 0; j < GRIDRES; j++)
193 –                if (fill_cnt[i][j])
194 –                        dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
195 – }
196 –
197 – /* Cull points for more uniform distribution, leave all nval 0 or 1 */
198 – static void
199 – cull_values(void)
200 – {
201 –        FVECT   ovec0, ovec1;
202 –        double  maxang, maxang2;
203 –        int     i, j, ii, jj, r;
204 –                                                /* simple greedy algorithm */
205 –        for (i = 0; i < GRIDRES; i++)
206 –            for (j = 0; j < GRIDRES; j++) {
207 –                if (!dsf_grid[i][j].nval)
208 –                        continue;
209 –                if (!dsf_grid[i][j].crad)
210 –                        continue;               /* shouldn't happen */
211 –                ovec_from_pos(ovec0, i, j);
212 –                maxang = 2.*R2ANG(dsf_grid[i][j].crad);
213 –                                                /* clamp near horizon */
214 –                if (maxang > output_orient*ovec0[2])
215 –                        maxang = output_orient*ovec0[2];
216 –                r = maxang*(2.*GRIDRES/M_PI) + 1;
217 –                maxang2 = maxang*maxang;
218 –                for (ii = i-r; ii <= i+r; ii++) {
219 –                    if (ii < 0) continue;
220 –                    if (ii >= GRIDRES) break;
221 –                    for (jj = j-r; jj <= j+r; jj++) {
222 –                        if ((ii == i) & (jj == j))
223 –                                continue;       /* don't get self-absorbed */
224 –                        if (jj < 0) continue;
225 –                        if (jj >= GRIDRES) break;
226 –                        if (!dsf_grid[ii][jj].nval)
227 –                                continue;
228 –                        ovec_from_pos(ovec1, ii, jj);
229 –                        if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
230 –                                continue;
231 –                                                /* absorb sum */
232 –                        dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
233 –                        dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
234 –                                                /* keep value, though */
235 –                        dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
236 –                        dsf_grid[ii][jj].nval = 0;
237 –                    }
238 –                }
239 –            }
240 –                                                /* final averaging pass */
241 –        for (i = 0; i < GRIDRES; i++)
242 –            for (j = 0; j < GRIDRES; j++)
243 –                if (dsf_grid[i][j].nval > 1) {
244 –                        dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
245 –                        dsf_grid[i][j].nval = 1;
246 –                }
247 – }
248 –
249 – /* Compute minimum BSDF from histogram and clear it */
250 – static void
90   comp_bsdf_min()
91   {
92 <        int     cnt;
93 <        int     i, target;
92 >        unsigned long   cnt, target;
93 >        int             i;
94  
95          cnt = 0;
96          for (i = HISTLEN; i--; )
# Line 265 | Line 104 | comp_bsdf_min()
104          for (i = 0; cnt <= target; i++)
105                  cnt += bsdf_hist[i];
106          bsdf_min = histval(i-1);
268 –        memset(bsdf_hist, 0, sizeof(bsdf_hist));
107   }
108  
109 < /* Find n nearest sub-sampled neighbors to the given grid position */
109 > /* Determine if the given region is empty of grid samples */
110   static int
111 < get_neighbors(int neigh[][2], int n, const int i, const int j)
111 > empty_region(int x0, int x1, int y0, int y1)
112   {
113 <        int     k = 0;
114 <        int     r;
115 <                                                /* search concentric squares */
116 <        for (r = 1; r < GRIDRES; r++) {
117 <                int     ii, jj;
118 <                for (ii = i-r; ii <= i+r; ii++) {
119 <                        int     jstep = 1;
120 <                        if (ii < 0) continue;
121 <                        if (ii >= GRIDRES) break;
122 <                        if ((i-r < ii) & (ii < i+r))
123 <                                jstep = r<<1;
124 <                        for (jj = j-r; jj <= j+r; jj += jstep) {
125 <                                if (jj < 0) continue;
126 <                                if (jj >= GRIDRES) break;
127 <                                if (dsf_grid[ii][jj].nval) {
128 <                                        neigh[k][0] = ii;
129 <                                        neigh[k][1] = jj;
130 <                                        if (++k >= n)
131 <                                                return(n);
132 <                                }
133 <                        }
113 >        int     x, y;
114 >
115 >        for (x = x0; x < x1; x++)
116 >            for (y = y0; y < y1; y++)
117 >                if (dsf_grid[x][y].sum.n)
118 >                        return(0);
119 >        return(1);
120 > }
121 >
122 > /* Determine if the given region is smooth enough to be a single lobe */
123 > static int
124 > smooth_region(int x0, int x1, int y0, int y1)
125 > {
126 >        RREAL   rMtx[3][3];
127 >        FVECT   xvec;
128 >        double  A, B, C, nvs, sqerr;
129 >        int     x, y, n;
130 >                                        /* compute planar regression */
131 >        memset(rMtx, 0, sizeof(rMtx));
132 >        memset(xvec, 0, sizeof(xvec));
133 >        for (x = x0; x < x1; x++)
134 >            for (y = y0; y < y1; y++)
135 >                if ((n = dsf_grid[x][y].sum.n) > 0) {
136 >                        double  z = dsf_grid[x][y].sum.v;
137 >                        rMtx[0][0] += x*x*(double)n;
138 >                        rMtx[0][1] += x*y*(double)n;
139 >                        rMtx[0][2] += x*(double)n;
140 >                        rMtx[1][1] += y*y*(double)n;
141 >                        rMtx[1][2] += y*(double)n;
142 >                        rMtx[2][2] += (double)n;
143 >                        xvec[0] += x*z;
144 >                        xvec[1] += y*z;
145 >                        xvec[2] += z;
146                  }
147 <        }
148 <        return(k);
147 >        rMtx[1][0] = rMtx[0][1];
148 >        rMtx[2][0] = rMtx[0][2];
149 >        rMtx[2][1] = rMtx[1][2];
150 >        nvs = rMtx[2][2];
151 >        if (SDinvXform(rMtx, rMtx) != SDEnone)
152 >                return(1);              /* colinear values */
153 >        A = DOT(rMtx[0], xvec);
154 >        B = DOT(rMtx[1], xvec);
155 >        C = DOT(rMtx[2], xvec);
156 >        sqerr = 0.0;                    /* compute mean squared error */
157 >        for (x = x0; x < x1; x++)
158 >            for (y = y0; y < y1; y++)
159 >                if ((n = dsf_grid[x][y].sum.n) > 0) {
160 >                        double  d = A*x + B*y + C - dsf_grid[x][y].sum.v/n;
161 >                        sqerr += n*d*d;
162 >                }
163 >        if (sqerr <= nvs*SMOOTH_MSE)    /* below absolute MSE threshold? */
164 >                return(1);
165 >                                        /* OR below relative MSE threshold? */
166 >        return(sqerr*nvs <= xvec[2]*xvec[2]*SMOOTH_MSER);
167   }
168  
169 < /* Adjust coded radius for the given grid position based on neighborhood */
169 > /* Create new lobe based on integrated samples in region */
170   static int
171 < adj_coded_radius(const int i, const int j)
171 > create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
172   {
173 <        const double    rad0 = R2ANG(dsf_grid[i][j].crad);
174 <        const double    minrad = MINRSCA * rad0;
175 <        double          currad = MAXRSCA * rad0;
176 <        int             neigh[NNEIGH][2];
177 <        int             n;
178 <        FVECT           our_dir;
173 >        double  vtot = 0.0;
174 >        int     nv = 0;
175 >        double  wxsum = 0.0, wysum = 0.0, wtsum = 0.0;
176 >        double  rad;
177 >        int     x, y;
178 >                                        /* compute average for region */
179 >        for (x = x0; x < x1; x++)
180 >            for (y = y0; y < y1; y++)
181 >                if (dsf_grid[x][y].sum.n) {
182 >                    const double        v = dsf_grid[x][y].sum.v;
183 >                    const int           n = dsf_grid[x][y].sum.n;
184  
185 <        ovec_from_pos(our_dir, i, j);
186 <        n = get_neighbors(neigh, NNEIGH, i, j);
187 <        while (n--) {
188 <                FVECT   their_dir;
189 <                double  max_ratio, rad_ok2;
190 <                                                /* check our value at neighbor */
191 <                ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
192 <                max_ratio = MAXFRAC * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
193 <                                / dsf_grid[i][j].vsum;
194 <                if (max_ratio >= 1)
195 <                        continue;
196 <                rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
197 <                if (rad_ok2 >= currad*currad)
325 <                        continue;               /* value fraction OK */
326 <                currad = sqrt(rad_ok2);         /* else reduce lobe radius */
327 <                if (currad <= minrad)           /* limit how small we'll go */
328 <                        return(ANG2R(minrad));
185 >                    if (v > 0) {
186 >                        double  wt = v / (double)n;
187 >                        wxsum += wt * x;
188 >                        wysum += wt * y;
189 >                        wtsum += wt;
190 >                    }
191 >                    vtot += v;
192 >                    nv += n;
193 >                }
194 >        if (!nv) {
195 >                fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
196 >                                progname);
197 >                exit(1);
198          }
199 <        return(ANG2R(currad));                  /* encode selected radius */
199 >        if (vtot <= 0)                  /* only create positive lobes */
200 >                return(0);
201 >                                        /* peak value based on integral */
202 >        vtot *= (x1-x0)*(y1-y0)*(2.*M_PI/GRIDRES/GRIDRES)/(double)nv;
203 >        rad = (RSCA/(double)GRIDRES)*(x1-x0);
204 >        rvp->peak =  vtot / ((2.*M_PI) * rad*rad);
205 >        rvp->crad = ANG2R(rad);         /* put peak at centroid */
206 >        rvp->gx = (int)(wxsum/wtsum + .5);
207 >        rvp->gy = (int)(wysum/wtsum + .5);
208 >        return(1);
209   }
210  
211 + /* Recursive function to build radial basis function representation */
212 + static int
213 + build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, int y0, int y1)
214 + {
215 +        int     xmid = (x0+x1)>>1;
216 +        int     ymid = (y0+y1)>>1;
217 +        int     branched[4];
218 +        int     nadded, nleaves;
219 +                                        /* need to make this a leaf? */
220 +        if (empty_region(x0, xmid, y0, ymid) ||
221 +                        empty_region(xmid, x1, y0, ymid) ||
222 +                        empty_region(x0, xmid, ymid, y1) ||
223 +                        empty_region(xmid, x1, ymid, y1))
224 +                return(0);
225 +                                        /* add children (branches+leaves) */
226 +        if ((branched[0] = build_rbfrep(arp, np, x0, xmid, y0, ymid)) < 0)
227 +                return(-1);
228 +        if ((branched[1] = build_rbfrep(arp, np, xmid, x1, y0, ymid)) < 0)
229 +                return(-1);
230 +        if ((branched[2] = build_rbfrep(arp, np, x0, xmid, ymid, y1)) < 0)
231 +                return(-1);
232 +        if ((branched[3] = build_rbfrep(arp, np, xmid, x1, ymid, y1)) < 0)
233 +                return(-1);
234 +        nadded = branched[0] + branched[1] + branched[2] + branched[3];
235 +        nleaves = !branched[0] + !branched[1] + !branched[2] + !branched[3];
236 +        if (!nleaves)                   /* nothing but branches? */
237 +                return(nadded);
238 +                                        /* combine 4 leaves into 1? */
239 +        if ((nleaves == 4) & (x1-x0 <= MAX_RAD) &&
240 +                        smooth_region(x0, x1, y0, y1))
241 +                return(0);
242 +                                        /* need more array space? */
243 +        if ((*np+nleaves-1)>>RBFALLOCB != (*np-1)>>RBFALLOCB) {
244 +                *arp = (RBFVAL *)realloc(*arp,
245 +                                sizeof(RBFVAL)*(*np+nleaves-1+(1<<RBFALLOCB)));
246 +                if (*arp == NULL)
247 +                        return(-1);
248 +        }
249 +                                        /* create lobes for leaves */
250 +        if (!branched[0] && create_lobe(*arp+*np, x0, xmid, y0, ymid)) {
251 +                ++(*np); ++nadded;
252 +        }
253 +        if (!branched[1] && create_lobe(*arp+*np, xmid, x1, y0, ymid)) {
254 +                ++(*np); ++nadded;
255 +        }
256 +        if (!branched[2] && create_lobe(*arp+*np, x0, xmid, ymid, y1)) {
257 +                ++(*np); ++nadded;
258 +        }
259 +        if (!branched[3] && create_lobe(*arp+*np, xmid, x1, ymid, y1)) {
260 +                ++(*np); ++nadded;
261 +        }
262 +        return(nadded);
263 + }
264 +
265   /* Count up filled nodes and build RBF representation from current grid */
266   RBFNODE *
267 < make_rbfrep(void)
267 > make_rbfrep()
268   {
337 –        long    cradsum = 0, ocradsum = 0;
338 –        int     niter = 16;
339 –        double  lastVar, thisVar = 100.;
340 –        int     nn;
269          RBFNODE *newnode;
270 <        RBFVAL  *itera;
271 <        int     i, j;
344 <
345 < #ifdef DEBUG
346 < {
347 <        int     maxcnt = 0, nempty = 0;
348 <        long    cntsum = 0;
349 <        for (i = 0; i < GRIDRES; i++)
350 <            for (j = 0; j < GRIDRES; j++)
351 <                if (!dsf_grid[i][j].nval) {
352 <                        ++nempty;
353 <                } else {
354 <                        if (dsf_grid[i][j].nval > maxcnt)
355 <                                maxcnt = dsf_grid[i][j].nval;
356 <                        cntsum += dsf_grid[i][j].nval;
357 <                }
358 <        fprintf(stderr, "Average, maximum bin count: %d, %d (%.1f%% empty)\n",
359 <                        (int)(cntsum/((GRIDRES*GRIDRES)-nempty)), maxcnt,
360 <                        100./(GRIDRES*GRIDRES)*nempty);
361 < }
362 < #endif
363 <                                /* compute RBF radii */
364 <        compute_radii();
365 <                                /* coagulate lobes */
366 <        cull_values();
367 <        nn = 0;                 /* count selected bins */
368 <        for (i = 0; i < GRIDRES; i++)
369 <            for (j = 0; j < GRIDRES; j++)
370 <                nn += dsf_grid[i][j].nval;
270 >        RBFVAL  *rbfarr;
271 >        int     nn;
272                                  /* compute minimum BSDF */
273          comp_bsdf_min();
274 <                                /* allocate RBF array */
275 <        newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
274 >                                /* create RBF node list */
275 >        rbfarr = NULL; nn = 0;
276 >        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0) {
277 >                if (nn)
278 >                        goto memerr;
279 >                fprintf(stderr,
280 >                        "%s: warning - skipping bad incidence (%.1f,%.1f)\n",
281 >                                progname, theta_in_deg, phi_in_deg);
282 >                return(NULL);
283 >        }
284 >                                /* (re)allocate RBF array */
285 >        newnode = (RBFNODE *)realloc(rbfarr,
286 >                        sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
287          if (newnode == NULL)
288                  goto memerr;
289 +                                /* copy computed lobes into RBF node */
290 +        memmove(newnode->rbfa, newnode, sizeof(RBFVAL)*nn);
291          newnode->ord = -1;
292          newnode->next = NULL;
293          newnode->ejl = NULL;
# Line 381 | Line 295 | make_rbfrep(void)
295          newnode->invec[0] = cos((M_PI/180.)*phi_in_deg)*newnode->invec[2];
296          newnode->invec[1] = sin((M_PI/180.)*phi_in_deg)*newnode->invec[2];
297          newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
298 <        newnode->vtotal = 0;
298 >        newnode->vtotal = .0;
299          newnode->nrbf = nn;
300 <        nn = 0;                 /* fill RBF array */
301 <        for (i = 0; i < GRIDRES; i++)
302 <            for (j = 0; j < GRIDRES; j++)
389 <                if (dsf_grid[i][j].nval) {
390 <                        newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
391 <                        ocradsum += dsf_grid[i][j].crad;
392 <                        cradsum +=
393 <                        newnode->rbfa[nn].crad = adj_coded_radius(i, j);
394 <                        newnode->rbfa[nn].gx = i;
395 <                        newnode->rbfa[nn].gy = j;
396 <                        ++nn;
397 <                }
300 >                                /* compute sum for normalization */
301 >        while (nn-- > 0)
302 >                newnode->vtotal += rbf_volume(&newnode->rbfa[nn]);
303   #ifdef DEBUG
304 <        fprintf(stderr,
400 <        "Average radius reduced from %.2f to %.2f degrees for %d lobes\n",
401 <                        180./M_PI*MAXRSCA*R2ANG(ocradsum/newnode->nrbf),
402 <                        180./M_PI*R2ANG(cradsum/newnode->nrbf), newnode->nrbf);
403 < #endif
404 <                                /* iterate to improve interpolation accuracy */
405 <        itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
406 <        if (itera == NULL)
407 <                goto memerr;
408 <        memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
409 <        do {
410 <                double  dsum = 0, dsum2 = 0;
411 <                nn = 0;
412 <                for (i = 0; i < GRIDRES; i++)
413 <                    for (j = 0; j < GRIDRES; j++)
414 <                        if (dsf_grid[i][j].nval) {
415 <                                FVECT   odir;
416 <                                double  corr;
417 <                                ovec_from_pos(odir, i, j);
418 <                                itera[nn++].peak *= corr =
419 <                                        dsf_grid[i][j].vsum /
420 <                                                eval_rbfrep(newnode, odir);
421 <                                dsum += 1. - corr;
422 <                                dsum2 += (1.-corr)*(1.-corr);
423 <                        }
424 <                memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
425 <                lastVar = thisVar;
426 <                thisVar = dsum2/(double)nn;
427 < #ifdef DEBUG
428 <                fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
429 <                                        100.*dsum/(double)nn,
430 <                                        100.*sqrt(thisVar));
431 < #endif
432 <        } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
433 <
434 <        free(itera);
435 <        nn = 0;                 /* compute sum for normalization */
436 <        while (nn < newnode->nrbf)
437 <                newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
438 < #ifdef DEBUG
304 >        fprintf(stderr, "Built RBF with %d lobes\n", newnode->nrbf);
305          fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
306                          get_theta180(newnode->invec), get_phi360(newnode->invec),
307                          newnode->vtotal);
308   #endif
309          insert_dsf(newnode);
444 –
310          return(newnode);
311   memerr:
312          fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);

Diff Legend

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