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.14 by greg, Mon Oct 21 22:03:03 2013 UTC vs.
Revision 2.25 by greg, Sun Mar 30 00:19:09 2014 UTC

# Line 29 | Line 29 | static const char RCSid[] = "$Id$";
29   #include "bsdfrep.h"
30  
31   #ifndef RSCA
32 < #define RSCA            2.2             /* radius scaling factor (empirical) */
32 > #define RSCA            2.0             /* radius scaling factor (empirical) */
33   #endif
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.07            /* acceptable relative MSE */
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 <                                /* our loaded grid for this incident angle */
44 >                                        /* loaded grid or comparison DSFs */
45   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
46  
47   /* Start new DSF input grid */
# Line 72 | 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;
77 <        else if (!isDSF)
75 >        if (!isDSF)
76                  val *= ovec[2];         /* convert from BSDF to DSF */
77  
78                                          /* update BSDF histogram */
# Line 83 | 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   /* Compute minimum BSDF from histogram (does not clear) */
89   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 116 | Line 114 | empty_region(int x0, int x1, int y0, int y1)
114  
115          for (x = x0; x < x1; x++)
116              for (y = y0; y < y1; y++)
117 <                if (dsf_grid[x][y].nval)
117 >                if (dsf_grid[x][y].sum.n)
118                          return(0);
119          return(1);
120   }
# Line 134 | Line 132 | smooth_region(int x0, int x1, int y0, int y1)
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].nval) > 0) {
136 <                        double  z = dsf_grid[x][y].vsum;
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;
# Line 147 | Line 145 | smooth_region(int x0, int x1, int y0, int y1)
145                          xvec[2] += z;
146                  }
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(0);
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].nval) > 0) {
160 <                        double  d = A*x + B*y + C - dsf_grid[x][y].vsum/n;
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? */
# Line 168 | Line 167 | smooth_region(int x0, int x1, int y0, int y1)
167   }
168  
169   /* Create new lobe based on integrated samples in region */
170 < static void
170 > static int
171   create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
172   {
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 <                vtot += dsf_grid[x][y].vsum;
182 <                nv += dsf_grid[x][y].nval;
183 <            }
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 >                    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 +        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);
206 <        rvp->gx = (x0+x1)>>1;
207 <        rvp->gy = (y0+y1)>>1;
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 */
# Line 234 | Line 247 | build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, in
247                          return(-1);
248          }
249                                          /* create lobes for leaves */
250 <        if (!branched[0])
251 <                create_lobe(*arp+(*np)++, x0, xmid, y0, ymid);
252 <        if (!branched[1])
253 <                create_lobe(*arp+(*np)++, xmid, x1, y0, ymid);
254 <        if (!branched[2])
255 <                create_lobe(*arp+(*np)++, x0, xmid, ymid, y1);
256 <        if (!branched[3])
257 <                create_lobe(*arp+(*np)++, xmid, x1, ymid, y1);
258 <        nadded += nleaves;
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  
# Line 257 | Line 273 | make_rbfrep()
273          comp_bsdf_min();
274                                  /* create RBF node list */
275          rbfarr = NULL; nn = 0;
276 <        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0)
277 <                goto memerr;
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));
# Line 285 | Line 307 | make_rbfrep()
307                          newnode->vtotal);
308   #endif
309          insert_dsf(newnode);
288 –
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