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.12 by greg, Mon Oct 21 18:33:15 2013 UTC vs.
Revision 2.27 by greg, Fri Jan 29 16:21:55 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.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 +                                        /* 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 37 | 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 (val <= 0)                   /* truncate to zero */
62 <                val = 0;
63 <        else if (!isDSF)
64 <                val *= ovec[2];         /* convert from BSDF to DSF */
65 <
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]] += val[0];
141 >                spec_grid[1][pos[0]][pos[1]] += val[2];
142 >        }
143   }
144  
145   /* Compute minimum BSDF from histogram (does not clear) */
146   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 102 | Line 171 | empty_region(int x0, int x1, int y0, int y1)
171  
172          for (x = x0; x < x1; x++)
173              for (y = y0; y < y1; y++)
174 <                if (dsf_grid[x][y].nval)
174 >                if (dsf_grid[x][y].sum.n)
175                          return(0);
176          return(1);
177   }
# Line 120 | Line 189 | smooth_region(int x0, int x1, int y0, int y1)
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].nval) > 0) {
193 <                        double  z = dsf_grid[x][y].vsum;
194 <                        rMtx[0][0] += n*x*x;
195 <                        rMtx[0][1] += n*x*y;
196 <                        rMtx[0][2] += n*x;
197 <                        rMtx[1][1] += n*y*y;
198 <                        rMtx[1][2] += n*y;
199 <                        rMtx[2][2] += n;
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(0);
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].nval) > 0) {
217 <                        double  d = A*x + B*y + C - dsf_grid[x][y].vsum/n;
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 <                                        /* below relative MSE threshold? */
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 void
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 <                vtot += dsf_grid[x][y].vsum;
240 <                nv += dsf_grid[x][y].nval;
241 <            }
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);
276 <        rvp->gx = (x0+x1)>>1;
277 <        rvp->gy = (y0+y1)>>1;
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 */
# Line 209 | Line 306 | build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, in
306          if (!nleaves)                   /* nothing but branches? */
307                  return(nadded);
308                                          /* combine 4 leaves into 1? */
309 <        if (nleaves == 4 && x1-x0 < MAX_RAD && smooth_region(x0, x1, y0, y1))
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) {
# Line 219 | Line 317 | build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, in
317                          return(-1);
318          }
319                                          /* create lobes for leaves */
320 <        if (!branched[0])
321 <                create_lobe(*arp+(*np)++, x0, xmid, y0, ymid);
322 <        if (!branched[1])
323 <                create_lobe(*arp+(*np)++, xmid, x1, y0, ymid);
324 <        if (!branched[2])
325 <                create_lobe(*arp+(*np)++, x0, xmid, ymid, y1);
326 <        if (!branched[3])
327 <                create_lobe(*arp+(*np)++, xmid, x1, ymid, y1);
328 <        nadded += nleaves;
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  
# Line 242 | Line 343 | make_rbfrep()
343          comp_bsdf_min();
344                                  /* create RBF node list */
345          rbfarr = NULL; nn = 0;
346 <        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0)
347 <                goto memerr;
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));
# Line 270 | Line 377 | make_rbfrep()
377                          newnode->vtotal);
378   #endif
379          insert_dsf(newnode);
273 –
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