ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfrbf.c
Revision: 2.11
Committed: Sat Oct 19 00:11:50 2013 UTC (10 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.10: +61 -34 lines
Log Message:
Fixed overflow and tweaked culling operations

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdfrbf.c,v 2.10 2013/10/18 02:49:30 greg Exp $";
3 #endif
4 /*
5 * Radial basis function representation for BSDF data.
6 *
7 * G. Ward
8 */
9
10 #define _USE_MATH_DEFINES
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <string.h>
14 #include <math.h>
15 #include "bsdfrep.h"
16
17 #ifndef MINRSCA
18 #define MINRSCA 1.0 /* minimum radius scaling factor */
19 #endif
20 #ifndef MAXRSCA
21 #define MAXRSCA 2.7 /* maximum radius scaling factor */
22 #endif
23 #ifndef VARTHRESH
24 #define VARTHRESH 0.0015 /* culling variance threshold */
25 #endif
26 #ifndef DIFFMAX2
27 #define DIFFMAX2 (16.*VARTHRESH) /* maximum ignored sample variance */
28 #endif
29 #ifndef MAXFRAC
30 #define MAXFRAC 0.5 /* maximum contribution to neighbor */
31 #endif
32 #ifndef NNEIGH
33 #define NNEIGH 10 /* number of neighbors to consider */
34 #endif
35 /* our loaded grid for this incident angle */
36 GRIDVAL dsf_grid[GRIDRES][GRIDRES];
37
38 /* Start new DSF input grid */
39 void
40 new_bsdf_data(double new_theta, double new_phi)
41 {
42 if (!new_input_direction(new_theta, new_phi))
43 exit(1);
44 memset(dsf_grid, 0, sizeof(dsf_grid));
45 }
46
47 /* Add BSDF data point */
48 void
49 add_bsdf_data(double theta_out, double phi_out, double val, int isDSF)
50 {
51 FVECT ovec;
52 int pos[2];
53
54 if (!output_orient) /* check output orientation */
55 output_orient = 1 - 2*(theta_out > 90.);
56 else if (output_orient > 0 ^ theta_out < 90.) {
57 fputs("Cannot handle output angles on both sides of surface\n",
58 stderr);
59 exit(1);
60 }
61 ovec[2] = sin((M_PI/180.)*theta_out);
62 ovec[0] = cos((M_PI/180.)*phi_out) * ovec[2];
63 ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
64 ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
65
66 if (val <= 0) /* truncate to zero */
67 val = 0;
68 else if (!isDSF)
69 val *= ovec[2]; /* convert from BSDF to DSF */
70
71 /* update BSDF histogram */
72 if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
73 ++bsdf_hist[histndx(val/ovec[2])];
74
75 pos_from_vec(pos, ovec);
76
77 dsf_grid[pos[0]][pos[1]].vsum += val;
78 dsf_grid[pos[0]][pos[1]].nval++;
79 }
80
81 /* Compute radii for non-empty bins */
82 /* (distance to furthest empty bin for which non-empty test bin is closest) */
83 static void
84 compute_radii(void)
85 {
86 const int cradmin = ANG2R(.5*M_PI/GRIDRES);
87 unsigned int fill_grid[GRIDRES][GRIDRES];
88 unsigned short fill_cnt[GRIDRES][GRIDRES];
89 FVECT ovec0, ovec1;
90 double ang2, lastang2;
91 int r, i, j, jn, ii, jj, inear, jnear;
92
93 r = GRIDRES/2; /* proceed in zig-zag */
94 for (i = 0; i < GRIDRES; i++)
95 for (jn = 0; jn < GRIDRES; jn++) {
96 j = (i&1) ? jn : GRIDRES-1-jn;
97 if (dsf_grid[i][j].nval) /* find empty grid pos. */
98 continue;
99 ovec_from_pos(ovec0, i, j);
100 inear = jnear = -1; /* find nearest non-empty */
101 lastang2 = M_PI*M_PI;
102 for (ii = i-r; ii <= i+r; ii++) {
103 if (ii < 0) continue;
104 if (ii >= GRIDRES) break;
105 for (jj = j-r; jj <= j+r; jj++) {
106 if (jj < 0) continue;
107 if (jj >= GRIDRES) break;
108 if (!dsf_grid[ii][jj].nval)
109 continue;
110 ovec_from_pos(ovec1, ii, jj);
111 ang2 = 2. - 2.*DOT(ovec0,ovec1);
112 if (ang2 >= lastang2)
113 continue;
114 lastang2 = ang2;
115 inear = ii; jnear = jj;
116 }
117 }
118 if (inear < 0) {
119 fprintf(stderr,
120 "%s: Could not find non-empty neighbor!\n",
121 progname);
122 exit(1);
123 }
124 ang2 = sqrt(lastang2);
125 r = ANG2R(ang2); /* record if > previous */
126 if (r > dsf_grid[inear][jnear].crad)
127 dsf_grid[inear][jnear].crad = r;
128 /* next search radius */
129 r = ang2*(2.*GRIDRES/M_PI) + 3;
130 }
131 for (i = 0; i < GRIDRES; i++) /* grow radii where uniform */
132 for (j = 0; j < GRIDRES; j++) {
133 double midmean = 0.0;
134 int nsum = 0;
135 if (!dsf_grid[i][j].nval)
136 continue;
137 r = 1; /* avg. immediate neighbors */
138 for (ii = i-r; ii <= i+r; ii++) {
139 if (ii < 0) continue;
140 if (ii >= GRIDRES) break;
141 for (jj = j-r; jj <= j+r; jj++) {
142 if (jj < 0) continue;
143 if (jj >= GRIDRES) break;
144 midmean += dsf_grid[ii][jj].vsum;
145 nsum += dsf_grid[ii][jj].nval;
146 }
147 }
148 midmean /= (double)nsum;
149 while (++r < GRIDRES) { /* attempt to grow perimeter */
150 double diff2sum = 0.0;
151 nsum = 0;
152 for (ii = i-r; ii <= i+r; ii++) {
153 int jstep = 1;
154 if (ii < 0) continue;
155 if (ii >= GRIDRES) break;
156 if ((i-r < ii) & (ii < i+r))
157 jstep = r<<1;
158 for (jj = j-r; jj <= j+r; jj += jstep) {
159 double d2;
160 if (jj < 0) continue;
161 if (jj >= GRIDRES) break;
162 if (!dsf_grid[ii][jj].nval)
163 continue;
164 d2 = midmean - dsf_grid[ii][jj].vsum /
165 (double)dsf_grid[ii][jj].nval;
166 d2 *= d2;
167 if (d2 > DIFFMAX2*midmean*midmean)
168 goto escape;
169 diff2sum += d2;
170 ++nsum;
171 }
172 }
173 if (diff2sum > VARTHRESH*midmean*midmean*(double)nsum)
174 break;
175 }
176 escape: --r;
177 r = ANG2R(r*(M_PI/MAXRSCA/GRIDRES));
178 if (r < cradmin)
179 r = cradmin;
180 if (dsf_grid[i][j].crad < r)
181 dsf_grid[i][j].crad = r;
182 }
183 /* blur radii over hemisphere */
184 memset(fill_grid, 0, sizeof(fill_grid));
185 memset(fill_cnt, 0, sizeof(fill_cnt));
186 for (i = 0; i < GRIDRES; i++)
187 for (j = 0; j < GRIDRES; j++) {
188 if (!dsf_grid[i][j].nval)
189 continue; /* not part of this */
190 r = R2ANG(dsf_grid[i][j].crad)*(2.*MAXRSCA*GRIDRES/M_PI);
191 for (ii = i-r; ii <= i+r; ii++) {
192 if (ii < 0) continue;
193 if (ii >= GRIDRES) break;
194 for (jj = j-r; jj <= j+r; jj++) {
195 if (jj < 0) continue;
196 if (jj >= GRIDRES) break;
197 if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
198 continue;
199 fill_grid[ii][jj] += dsf_grid[i][j].crad;
200 fill_cnt[ii][jj]++;
201 }
202 }
203 }
204 /* copy back blurred radii */
205 for (i = 0; i < GRIDRES; i++)
206 for (j = 0; j < GRIDRES; j++)
207 if (fill_cnt[i][j])
208 dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
209 }
210
211 /* Radius comparison for qsort() */
212 static int
213 radius_cmp(const void *p1, const void *p2)
214 {
215 return( (int)dsf_grid[0][*(const int *)p1].crad -
216 (int)dsf_grid[0][*(const int *)p2].crad );
217 }
218
219 /* Cull points for more uniform distribution, leave all nval 0 or 1 */
220 static void
221 cull_values(void)
222 {
223 int indx[GRIDRES*GRIDRES];
224 FVECT ovec0, ovec1;
225 double maxang, maxang2;
226 int i, j, k, ii, jj, r;
227 /* sort by radius first */
228 for (k = GRIDRES*GRIDRES; k--; )
229 indx[k] = k;
230 qsort(indx, GRIDRES*GRIDRES, sizeof(int), &radius_cmp);
231 /* simple greedy algorithm */
232 for (k = GRIDRES*GRIDRES; k--; ) {
233 i = indx[k]/GRIDRES; /* from biggest radius down */
234 j = indx[k] - i*GRIDRES;
235 if (!dsf_grid[i][j].nval)
236 continue;
237 if (!dsf_grid[i][j].crad)
238 break; /* shouldn't happen */
239 ovec_from_pos(ovec0, i, j);
240 maxang = 2.*R2ANG(dsf_grid[i][j].crad);
241 /* clamp near horizon */
242 if (maxang > output_orient*ovec0[2])
243 maxang = output_orient*ovec0[2];
244 r = maxang*(2.*GRIDRES/M_PI) + 1;
245 maxang2 = maxang*maxang;
246 for (ii = i-r; ii <= i+r; ii++) {
247 if (ii < 0) continue;
248 if (ii >= GRIDRES) break;
249 for (jj = j-r; jj <= j+r; jj++) {
250 if ((ii == i) & (jj == j))
251 continue; /* don't get self-absorbed */
252 if (jj < 0) continue;
253 if (jj >= GRIDRES) break;
254 if (!dsf_grid[ii][jj].nval)
255 continue;
256 ovec_from_pos(ovec1, ii, jj);
257 if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
258 continue;
259 /* absorb sum */
260 dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
261 dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
262 /* keep value, though */
263 dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
264 dsf_grid[ii][jj].nval = 0;
265 }
266 }
267 }
268 /* final averaging pass */
269 for (i = 0; i < GRIDRES; i++)
270 for (j = 0; j < GRIDRES; j++)
271 if (dsf_grid[i][j].nval > 1) {
272 dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
273 dsf_grid[i][j].nval = 1;
274 }
275 }
276
277 /* Compute minimum BSDF from histogram (does not clear) */
278 static void
279 comp_bsdf_min()
280 {
281 int cnt;
282 int i, target;
283
284 cnt = 0;
285 for (i = HISTLEN; i--; )
286 cnt += bsdf_hist[i];
287 if (!cnt) { /* shouldn't happen */
288 bsdf_min = 0;
289 return;
290 }
291 target = cnt/100; /* ignore bottom 1% */
292 cnt = 0;
293 for (i = 0; cnt <= target; i++)
294 cnt += bsdf_hist[i];
295 bsdf_min = histval(i-1);
296 }
297
298 /* Find n nearest sub-sampled neighbors to the given grid position */
299 static int
300 get_neighbors(int neigh[][2], int n, const int i, const int j)
301 {
302 int k = 0;
303 int r;
304 /* search concentric squares */
305 for (r = 1; r < GRIDRES; r++) {
306 int ii, jj;
307 for (ii = i-r; ii <= i+r; ii++) {
308 int jstep = 1;
309 if (ii < 0) continue;
310 if (ii >= GRIDRES) break;
311 if ((i-r < ii) & (ii < i+r))
312 jstep = r<<1;
313 for (jj = j-r; jj <= j+r; jj += jstep) {
314 if (jj < 0) continue;
315 if (jj >= GRIDRES) break;
316 if (dsf_grid[ii][jj].nval) {
317 neigh[k][0] = ii;
318 neigh[k][1] = jj;
319 if (++k >= n)
320 return(n);
321 }
322 }
323 }
324 }
325 return(k);
326 }
327
328 /* Adjust coded radius for the given grid position based on neighborhood */
329 static int
330 adj_coded_radius(const int i, const int j)
331 {
332 const double rad0 = R2ANG(dsf_grid[i][j].crad);
333 const double minrad = MINRSCA * rad0;
334 double currad = MAXRSCA * rad0;
335 int neigh[NNEIGH][2];
336 int n;
337 FVECT our_dir;
338
339 ovec_from_pos(our_dir, i, j);
340 n = get_neighbors(neigh, NNEIGH, i, j);
341 while (n--) {
342 FVECT their_dir;
343 double max_ratio, rad_ok2;
344 /* check our value at neighbor */
345 ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
346 max_ratio = MAXFRAC * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
347 / dsf_grid[i][j].vsum;
348 if (max_ratio >= 1)
349 continue;
350 rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
351 if (rad_ok2 >= currad*currad)
352 continue; /* value fraction OK */
353 currad = sqrt(rad_ok2); /* else reduce lobe radius */
354 if (currad <= minrad) /* limit how small we'll go */
355 return(ANG2R(minrad));
356 }
357 return(ANG2R(currad)); /* encode selected radius */
358 }
359
360 /* Count up filled nodes and build RBF representation from current grid */
361 RBFNODE *
362 make_rbfrep(void)
363 {
364 long cradsum = 0, ocradsum = 0;
365 int niter = 16;
366 double lastVar, thisVar = 100.;
367 int nn;
368 RBFNODE *newnode;
369 RBFVAL *itera;
370 int i, j;
371
372 #ifdef DEBUG
373 {
374 int maxcnt = 0, nempty = 0;
375 long cntsum = 0;
376 for (i = 0; i < GRIDRES; i++)
377 for (j = 0; j < GRIDRES; j++)
378 if (!dsf_grid[i][j].nval) {
379 ++nempty;
380 } else {
381 if (dsf_grid[i][j].nval > maxcnt)
382 maxcnt = dsf_grid[i][j].nval;
383 cntsum += dsf_grid[i][j].nval;
384 }
385 fprintf(stderr, "Average, maximum bin count: %d, %d (%.1f%% empty)\n",
386 (int)(cntsum/((GRIDRES*GRIDRES)-nempty)), maxcnt,
387 100./(GRIDRES*GRIDRES)*nempty);
388 }
389 #endif
390 /* compute RBF radii */
391 compute_radii();
392 /* coagulate lobes */
393 cull_values();
394 nn = 0; /* count selected bins */
395 for (i = 0; i < GRIDRES; i++)
396 for (j = 0; j < GRIDRES; j++)
397 nn += dsf_grid[i][j].nval;
398 /* compute minimum BSDF */
399 comp_bsdf_min();
400 /* allocate RBF array */
401 newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
402 if (newnode == NULL)
403 goto memerr;
404 newnode->ord = -1;
405 newnode->next = NULL;
406 newnode->ejl = NULL;
407 newnode->invec[2] = sin((M_PI/180.)*theta_in_deg);
408 newnode->invec[0] = cos((M_PI/180.)*phi_in_deg)*newnode->invec[2];
409 newnode->invec[1] = sin((M_PI/180.)*phi_in_deg)*newnode->invec[2];
410 newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
411 newnode->vtotal = 0;
412 newnode->nrbf = nn;
413 nn = 0; /* fill RBF array */
414 for (i = 0; i < GRIDRES; i++)
415 for (j = 0; j < GRIDRES; j++)
416 if (dsf_grid[i][j].nval) {
417 newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
418 ocradsum += dsf_grid[i][j].crad;
419 cradsum +=
420 newnode->rbfa[nn].crad = adj_coded_radius(i, j);
421 newnode->rbfa[nn].gx = i;
422 newnode->rbfa[nn].gy = j;
423 ++nn;
424 }
425 #ifdef DEBUG
426 fprintf(stderr,
427 "Average radius reduced from %.2f to %.2f degrees for %d lobes\n",
428 180./M_PI*MAXRSCA*R2ANG(ocradsum/newnode->nrbf),
429 180./M_PI*R2ANG(cradsum/newnode->nrbf), newnode->nrbf);
430 #endif
431 /* iterate to improve interpolation accuracy */
432 itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
433 if (itera == NULL)
434 goto memerr;
435 memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
436 do {
437 double dsum = 0, dsum2 = 0;
438 nn = 0;
439 for (i = 0; i < GRIDRES; i++)
440 for (j = 0; j < GRIDRES; j++)
441 if (dsf_grid[i][j].nval) {
442 FVECT odir;
443 double corr;
444 ovec_from_pos(odir, i, j);
445 itera[nn++].peak *= corr =
446 dsf_grid[i][j].vsum /
447 eval_rbfrep(newnode, odir);
448 dsum += 1. - corr;
449 dsum2 += (1.-corr)*(1.-corr);
450 }
451 memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
452 lastVar = thisVar;
453 thisVar = dsum2/(double)nn;
454 #ifdef DEBUG
455 fprintf(stderr, "Avg., RMS error: %.1f%% %.1f%%\n",
456 100.*dsum/(double)nn,
457 100.*sqrt(thisVar));
458 #endif
459 } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
460
461 free(itera);
462 nn = 0; /* compute sum for normalization */
463 while (nn < newnode->nrbf)
464 newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
465 #ifdef DEBUG
466 fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
467 get_theta180(newnode->invec), get_phi360(newnode->invec),
468 newnode->vtotal);
469 #endif
470 insert_dsf(newnode);
471
472 return(newnode);
473 memerr:
474 fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
475 exit(1);
476 }