ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfrep.c
Revision: 2.38
Committed: Sun May 18 01:46:05 2025 UTC (10 days, 17 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.37: +47 -3 lines
Log Message:
perf(bsdf2klems,bsdf2ttree): Improved evaluation speed by 40% using 128K table

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdfrep.c,v 2.37 2021/12/15 01:38:50 greg Exp $";
3 #endif
4 /*
5 * Support BSDF representation as radial basis functions.
6 *
7 * G. Ward
8 */
9
10 #define _USE_MATH_DEFINES
11 #include <stdlib.h>
12 #include <math.h>
13 #include "rtio.h"
14 #include "resolu.h"
15 #include "bsdfrep.h"
16 #include "random.h"
17 /* name and manufacturer if known */
18 char bsdf_name[256];
19 char bsdf_manuf[256];
20 /* active grid resolution */
21 int grid_res = GRIDRES;
22
23 /* coverage/symmetry using INP_QUAD? flags */
24 int inp_coverage = 0;
25 /* all incident angles in-plane so far? */
26 int single_plane_incident = -1;
27
28 /* input/output orientations */
29 int input_orient = 0;
30 int output_orient = 0;
31
32 /* represented color space */
33 RBColor rbf_colorimetry = RBCunknown;
34
35 const char *RBCident[] = {
36 "CIE-Y", "CIE-XYZ", "Spectral", "Unknown"
37 };
38
39 /* BSDF histogram */
40 unsigned long bsdf_hist[HISTLEN];
41
42 /* BSDF value for boundary regions */
43 double bsdf_min = 0;
44 double bsdf_spec_val = 0;
45 double bsdf_spec_rad = 0;
46
47 /* processed incident DSF measurements */
48 RBFNODE *dsf_list = NULL;
49
50 /* RBF-linking matrices (edges) */
51 MIGRATION *mig_list = NULL;
52
53 /* current input direction */
54 double theta_in_deg, phi_in_deg;
55
56 /* header line sharing callback */
57 int (*sir_headshare)(char *s) = NULL;
58
59 /* Register new input direction */
60 int
61 new_input_direction(double new_theta, double new_phi)
62 {
63 /* normalize angle ranges */
64 while (new_theta < -180.)
65 new_theta += 360.;
66 while (new_theta > 180.)
67 new_theta -= 360.;
68 if (new_theta < 0) {
69 new_theta = -new_theta;
70 new_phi += 180.;
71 }
72 while (new_phi < 0)
73 new_phi += 360.;
74 while (new_phi >= 360.)
75 new_phi -= 360.;
76 /* check input orientation */
77 if (!input_orient)
78 input_orient = 1 - 2*(new_theta > 90.);
79 else if (input_orient > 0 ^ new_theta < 90.) {
80 fprintf(stderr,
81 "%s: Cannot handle input angles on both sides of surface\n",
82 progname);
83 return(0);
84 }
85 if ((theta_in_deg = new_theta) < 1.0)
86 return(1); /* don't rely on phi near normal */
87 if (single_plane_incident > 0) /* check input coverage */
88 single_plane_incident = (round(new_phi) == round(phi_in_deg));
89 else if (single_plane_incident < 0)
90 single_plane_incident = 1;
91 phi_in_deg = new_phi;
92 if ((1. < new_phi) & (new_phi < 89.))
93 inp_coverage |= INP_QUAD1;
94 else if ((91. < new_phi) & (new_phi < 179.))
95 inp_coverage |= INP_QUAD2;
96 else if ((181. < new_phi) & (new_phi < 269.))
97 inp_coverage |= INP_QUAD3;
98 else if ((271. < new_phi) & (new_phi < 359.))
99 inp_coverage |= INP_QUAD4;
100 return(1);
101 }
102
103 /* Apply symmetry to the given vector based on distribution */
104 int
105 use_symmetry(FVECT vec)
106 {
107 double phi = get_phi360(vec);
108 /* because of -0. issue */
109 while (phi >= 360.) phi -= 360.;
110 while (phi < 0.) phi += 360.;
111
112 switch (inp_coverage) {
113 case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
114 break;
115 case INP_QUAD1|INP_QUAD2:
116 if ((-FTINY > phi) | (phi > 180.+FTINY))
117 goto mir_y;
118 break;
119 case INP_QUAD2|INP_QUAD3:
120 if ((90.-FTINY > phi) | (phi > 270.+FTINY))
121 goto mir_x;
122 break;
123 case INP_QUAD3|INP_QUAD4:
124 if ((180.-FTINY > phi) | (phi > 360.+FTINY))
125 goto mir_y;
126 break;
127 case INP_QUAD4|INP_QUAD1:
128 if ((270.-FTINY > phi) & (phi > 90.+FTINY))
129 goto mir_x;
130 break;
131 case INP_QUAD1:
132 if ((-FTINY > phi) | (phi > 90.+FTINY))
133 switch ((int)(phi*(1./90.))) {
134 case 1: goto mir_x;
135 case 2: goto mir_xy;
136 case 3: goto mir_y;
137 }
138 break;
139 case INP_QUAD2:
140 if ((90.-FTINY > phi) | (phi > 180.+FTINY))
141 switch ((int)(phi*(1./90.))) {
142 case 0: goto mir_x;
143 case 2: goto mir_y;
144 case 3: goto mir_xy;
145 }
146 break;
147 case INP_QUAD3:
148 if ((180.-FTINY > phi) | (phi > 270.+FTINY))
149 switch ((int)(phi*(1./90.))) {
150 case 0: goto mir_xy;
151 case 1: goto mir_y;
152 case 3: goto mir_x;
153 }
154 break;
155 case INP_QUAD4:
156 if ((270.-FTINY > phi) | (phi > 360.+FTINY))
157 switch ((int)(phi*(1./90.))) {
158 case 0: goto mir_y;
159 case 1: goto mir_xy;
160 case 2: goto mir_x;
161 }
162 break;
163 default:
164 fprintf(stderr, "%s: Illegal input coverage (%d)\n",
165 progname, inp_coverage);
166 exit(1);
167 }
168 return(0); /* in range */
169 mir_x:
170 vec[0] = -vec[0];
171 return(MIRROR_X);
172 mir_y:
173 vec[1] = -vec[1];
174 return(MIRROR_Y);
175 mir_xy:
176 vec[0] = -vec[0];
177 vec[1] = -vec[1];
178 return(MIRROR_X|MIRROR_Y);
179 }
180
181 /* Reverse symmetry based on what was done before */
182 void
183 rev_symmetry(FVECT vec, int sym)
184 {
185 if (sym & MIRROR_X)
186 vec[0] = -vec[0];
187 if (sym & MIRROR_Y)
188 vec[1] = -vec[1];
189 }
190
191 /* Reverse symmetry for an RBF distribution */
192 void
193 rev_rbf_symmetry(RBFNODE *rbf, int sym)
194 {
195 int n;
196
197 rev_symmetry(rbf->invec, sym);
198 if (sym & MIRROR_X)
199 for (n = rbf->nrbf; n-- > 0; )
200 rbf->rbfa[n].gx = grid_res-1 - rbf->rbfa[n].gx;
201 if (sym & MIRROR_Y)
202 for (n = rbf->nrbf; n-- > 0; )
203 rbf->rbfa[n].gy = grid_res-1 - rbf->rbfa[n].gy;
204 }
205
206 /* Rotate RBF to correspond to given incident vector */
207 void
208 rotate_rbf(RBFNODE *rbf, const FVECT invec)
209 {
210 static const FVECT vnorm = {.0, .0, 1.};
211 const double phi = atan2(invec[1],invec[0]) -
212 atan2(rbf->invec[1],rbf->invec[0]);
213 FVECT outvec;
214 int pos[2];
215 int n;
216
217 for (n = (cos(phi) < 1.-FTINY)*rbf->nrbf; n-- > 0; ) {
218 ovec_from_pos(outvec, rbf->rbfa[n].gx, rbf->rbfa[n].gy);
219 spinvector(outvec, outvec, vnorm, phi);
220 pos_from_vec(pos, outvec);
221 rbf->rbfa[n].gx = pos[0];
222 rbf->rbfa[n].gy = pos[1];
223 }
224 VCOPY(rbf->invec, invec);
225 }
226
227 /* Compute outgoing vector from grid position */
228 #if 1
229 void
230 ovec_from_pos(FVECT vec, int xpos, int ypos)
231 { /* precomputed table version */
232 static int qsiz = 0;
233 static float (*q_uv)[2] = NULL;
234
235 if (vec == NULL) { /* just free table? */
236 if (q_uv) free(q_uv);
237 qsiz = 0;
238 return;
239 }
240 if (qsiz != grid_res>>1) {
241 int x, y; /* (re)make positive quadrant table */
242 RREAL uv[2];
243 double r;
244 if (q_uv) free(q_uv);
245 qsiz = grid_res>>1;
246 q_uv = (float (*)[2])malloc(sizeof(float)*2*qsiz*qsiz);
247 for (y = qsiz; y--; )
248 for (x = qsiz; x--; ) {
249 square2disk(uv, 0.5 + (x+.5)/grid_res,
250 0.5 + (y+.5)/grid_res);
251 /* uniform hemispherical projection */
252 r = sqrt(2. - uv[0]*uv[0] - uv[1]*uv[1]);
253 q_uv[qsiz*y + x][0] = (float)(r*uv[0]);
254 q_uv[qsiz*y + x][1] = (float)(r*uv[1]);
255 }
256 }
257 /* put in positive quadrant */
258 if (xpos >= qsiz) { xpos -= qsiz; vec[0] = 1.; }
259 else { xpos = qsiz-1 - xpos; vec[0] = -1.; }
260 if (ypos >= qsiz) { ypos -= qsiz; vec[1] = 1.; }
261 else { ypos = qsiz-1 - ypos; vec[1] = -1.; }
262
263 vec[0] *= (RREAL)q_uv[qsiz*ypos + xpos][0];
264 vec[1] *= (RREAL)q_uv[qsiz*ypos + xpos][1];
265 vec[2] = output_orient*sqrt(1. - vec[0]*vec[0] - vec[1]*vec[1]);
266 }
267 #else
268 void
269 ovec_from_pos(FVECT vec, int xpos, int ypos)
270 { /* table-free version */
271 RREAL uv[2];
272 double r2;
273
274 if (vec == NULL)
275 return;
276
277 square2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
278 /* uniform hemispherical projection */
279 r2 = uv[0]*uv[0] + uv[1]*uv[1];
280 vec[0] = vec[1] = sqrt(2. - r2);
281 vec[0] *= uv[0];
282 vec[1] *= uv[1];
283 vec[2] = output_orient*(1. - r2);
284 }
285 #endif
286
287 /* Compute grid position from normalized input/output vector */
288 void
289 pos_from_vec(int pos[2], const FVECT vec)
290 {
291 RREAL sq[2]; /* uniform hemispherical projection */
292 double norm = 1./sqrt(1. + fabs(vec[2]));
293
294 disk2square(sq, vec[0]*norm, vec[1]*norm);
295
296 pos[0] = (int)(sq[0]*grid_res);
297 pos[1] = (int)(sq[1]*grid_res);
298 }
299
300 /* Compute volume associated with Gaussian lobe */
301 double
302 rbf_volume(const RBFVAL *rbfp)
303 {
304 double rad = R2ANG(rbfp->crad);
305 FVECT odir;
306 double elev, integ;
307 /* infinite integral approximation */
308 integ = (2.*M_PI) * rbfp->peak * rad*rad;
309 /* check if we're near horizon */
310 ovec_from_pos(odir, rbfp->gx, rbfp->gy);
311 elev = output_orient*odir[2];
312 /* apply cut-off correction if > 1% */
313 if (elev < 2.8*rad) {
314 /* elev = asin(elev); /* this is so crude, anyway... */
315 integ *= 1. - .5*exp(-.5*elev*elev/(rad*rad));
316 }
317 return(integ);
318 }
319
320 /* Evaluate BSDF at the given normalized outgoing direction in color */
321 SDError
322 eval_rbfcol(SDValue *sv, const RBFNODE *rp, const FVECT outvec)
323 {
324 const double rfact2 = (38./M_PI/M_PI)*(grid_res*grid_res);
325 int pos[2];
326 double res = 0;
327 double usum = 0, vsum = 0;
328 const RBFVAL *rbfp;
329 FVECT odir;
330 double rad2;
331 int n;
332 /* assign default value */
333 sv->spec = c_dfcolor;
334 sv->cieY = bsdf_min;
335 /* check for wrong side */
336 if (outvec[2] > 0 ^ output_orient > 0) {
337 strcpy(SDerrorDetail, "Wrong-side scattering query");
338 return(SDEargument);
339 }
340 if (rp == NULL) /* return minimum if no information avail. */
341 return(SDEnone);
342 /* optimization for fast lobe culling */
343 pos_from_vec(pos, outvec);
344 /* sum radial basis function */
345 rbfp = rp->rbfa;
346 for (n = rp->nrbf; n--; rbfp++) {
347 int d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
348 (pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
349 double val;
350 rad2 = R2ANG(rbfp->crad);
351 rad2 *= rad2;
352 if (d2 > rad2*rfact2)
353 continue;
354 ovec_from_pos(odir, rbfp->gx, rbfp->gy);
355 val = rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
356 if (rbf_colorimetry == RBCtristimulus) {
357 usum += val * (rbfp->chroma & 0xff);
358 vsum += val * (rbfp->chroma>>8 & 0xff);
359 }
360 res += val;
361 }
362 sv->cieY = res / COSF(outvec[2]);
363 if (sv->cieY < bsdf_min) { /* never return less than bsdf_min */
364 sv->cieY = bsdf_min;
365 } else if (rbf_colorimetry == RBCtristimulus) {
366 C_CHROMA cres = (int)(usum/res + frandom());
367 cres |= (int)(vsum/res + frandom()) << 8;
368 c_decodeChroma(&sv->spec, cres);
369 }
370 return(SDEnone);
371 }
372
373 /* Evaluate BSDF at the given normalized outgoing direction in Y */
374 double
375 eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
376 {
377 SDValue sv;
378
379 if (eval_rbfcol(&sv, rp, outvec) == SDEnone)
380 return(sv.cieY);
381
382 return(0.0);
383 }
384
385 /* Insert a new directional scattering function in our global list */
386 int
387 insert_dsf(RBFNODE *newrbf)
388 {
389 RBFNODE *rbf, *rbf_last;
390 int pos;
391 /* check for redundant meas. */
392 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
393 if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
394 fprintf(stderr,
395 "%s: Duplicate incident measurement ignored at (%.1f,%.1f)\n",
396 progname, get_theta180(newrbf->invec),
397 get_phi360(newrbf->invec));
398 free(newrbf);
399 return(-1);
400 }
401 /* keep in ascending theta order */
402 for (rbf_last = NULL, rbf = dsf_list; rbf != NULL;
403 rbf_last = rbf, rbf = rbf->next)
404 if (single_plane_incident && input_orient*rbf->invec[2] <
405 input_orient*newrbf->invec[2])
406 break;
407 if (rbf_last == NULL) { /* insert new node in list */
408 newrbf->ord = 0;
409 newrbf->next = dsf_list;
410 dsf_list = newrbf;
411 } else {
412 newrbf->ord = rbf_last->ord + 1;
413 newrbf->next = rbf;
414 rbf_last->next = newrbf;
415 }
416 rbf_last = newrbf;
417 while (rbf != NULL) { /* update ordinal positions */
418 rbf->ord = rbf_last->ord + 1;
419 rbf_last = rbf;
420 rbf = rbf->next;
421 }
422 return(newrbf->ord);
423 }
424
425 /* Get the DSF indicated by its ordinal position */
426 RBFNODE *
427 get_dsf(int ord)
428 {
429 RBFNODE *rbf;
430
431 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
432 if (rbf->ord == ord)
433 return(rbf);
434 return(NULL);
435 }
436
437 /* Get triangle surface orientation (unnormalized) */
438 void
439 tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
440 {
441 FVECT v2minus1, v3minus2;
442
443 VSUB(v2minus1, v2, v1);
444 VSUB(v3minus2, v3, v2);
445 VCROSS(vres, v2minus1, v3minus2);
446 }
447
448 /* Determine if vertex order is reversed (inward normal) */
449 int
450 is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
451 {
452 FVECT tor;
453
454 tri_orient(tor, v1, v2, v3);
455
456 return(DOT(tor, v2) < 0.);
457 }
458
459 /* Find vertices completing triangles on either side of the given edge */
460 int
461 get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
462 {
463 const MIGRATION *ej1, *ej2;
464 RBFNODE *tv;
465
466 rbfv[0] = rbfv[1] = NULL;
467 if (mig == NULL)
468 return(0);
469 for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
470 ej1 = nextedge(mig->rbfv[0],ej1)) {
471 if (ej1 == mig)
472 continue;
473 tv = opp_rbf(mig->rbfv[0],ej1);
474 for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
475 if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
476 rbfv[is_rev_tri(mig->rbfv[0]->invec,
477 mig->rbfv[1]->invec,
478 tv->invec)] = tv;
479 break;
480 }
481 }
482 return((rbfv[0] != NULL) + (rbfv[1] != NULL));
483 }
484
485 /* Return single-lobe specular RBF for the given incident direction */
486 RBFNODE *
487 def_rbf_spec(const FVECT invec)
488 {
489 RBFNODE *rbf;
490 FVECT ovec;
491 int pos[2];
492
493 if (input_orient > 0 ^ invec[2] > 0) /* wrong side? */
494 return(NULL);
495 if ((bsdf_spec_val <= bsdf_min) | (bsdf_spec_rad <= 0))
496 return(NULL); /* nothing set */
497 rbf = (RBFNODE *)malloc(sizeof(RBFNODE));
498 if (rbf == NULL)
499 return(NULL);
500 ovec[0] = -invec[0];
501 ovec[1] = -invec[1];
502 ovec[2] = invec[2]*(2*(input_orient==output_orient) - 1);
503 pos_from_vec(pos, ovec);
504 rbf->ord = 0;
505 rbf->next = NULL;
506 rbf->ejl = NULL;
507 VCOPY(rbf->invec, invec);
508 rbf->nrbf = 1;
509 rbf->rbfa[0].peak = bsdf_spec_val * COSF(ovec[2]);
510 rbf->rbfa[0].chroma = c_dfchroma;
511 rbf->rbfa[0].crad = ANG2R(bsdf_spec_rad);
512 rbf->rbfa[0].gx = pos[0];
513 rbf->rbfa[0].gy = pos[1];
514 rbf->vtotal = rbf_volume(rbf->rbfa);
515 return(rbf);
516 }
517
518 /* Advect and allocate new RBF along edge (internal call) */
519 RBFNODE *
520 e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
521 {
522 double cthresh = FTINY;
523 RBFNODE *rbf;
524 int n, i, j;
525 double t, full_dist;
526 /* get relative position */
527 t = Acos(DOT(invec, mig->rbfv[0]->invec));
528 if (t <= .001) { /* near first DSF */
529 n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
530 rbf = (RBFNODE *)malloc(n);
531 if (rbf == NULL)
532 goto memerr;
533 memcpy(rbf, mig->rbfv[0], n); /* just duplicate */
534 rbf->next = NULL; rbf->ejl = NULL;
535 return(rbf);
536 }
537 full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
538 if (t >= full_dist-.001) { /* near second DSF */
539 n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
540 rbf = (RBFNODE *)malloc(n);
541 if (rbf == NULL)
542 goto memerr;
543 memcpy(rbf, mig->rbfv[1], n); /* just duplicate */
544 rbf->next = NULL; rbf->ejl = NULL;
545 return(rbf);
546 }
547 t /= full_dist;
548 tryagain:
549 n = 0; /* count migrating particles */
550 for (i = 0; i < mtx_nrows(mig); i++)
551 for (j = 0; j < mtx_ncols(mig); j++)
552 n += (mtx_coef(mig,i,j) > cthresh);
553 /* are we over our limit? */
554 if ((lobe_lim > 0) & (n > lobe_lim)) {
555 cthresh = cthresh*2. + 10.*FTINY;
556 goto tryagain;
557 }
558 #ifdef DEBUG
559 fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
560 mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
561 #endif
562 rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
563 if (rbf == NULL)
564 goto memerr;
565 rbf->next = NULL; rbf->ejl = NULL;
566 VCOPY(rbf->invec, invec);
567 rbf->nrbf = n;
568 rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
569 n = 0; /* advect RBF lobes */
570 for (i = 0; i < mtx_nrows(mig); i++) {
571 const RBFVAL *rbf0i = &mig->rbfv[0]->rbfa[i];
572 const float peak0 = rbf0i->peak;
573 const double rad0 = R2ANG(rbf0i->crad);
574 C_COLOR cc0;
575 FVECT v0;
576 float mv;
577 ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
578 c_decodeChroma(&cc0, rbf0i->chroma);
579 for (j = 0; j < mtx_ncols(mig); j++)
580 if ((mv = mtx_coef(mig,i,j)) > cthresh) {
581 const RBFVAL *rbf1j = &mig->rbfv[1]->rbfa[j];
582 double rad2;
583 FVECT v;
584 int pos[2];
585 rad2 = R2ANG(rbf1j->crad);
586 rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
587 rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
588 rad0*rad0/rad2;
589 if (rbf_colorimetry == RBCtristimulus) {
590 C_COLOR cres;
591 c_decodeChroma(&cres, rbf1j->chroma);
592 c_cmix(&cres, 1.-t, &cc0, t, &cres);
593 rbf->rbfa[n].chroma = c_encodeChroma(&cres);
594 } else
595 rbf->rbfa[n].chroma = c_dfchroma;
596 rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
597 ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
598 geodesic(v, v0, v, t, GEOD_REL);
599 pos_from_vec(pos, v);
600 rbf->rbfa[n].gx = pos[0];
601 rbf->rbfa[n].gy = pos[1];
602 ++n;
603 }
604 }
605 rbf->vtotal *= mig->rbfv[0]->vtotal; /* turn ratio into actual */
606 return(rbf);
607 memerr:
608 fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname);
609 exit(1);
610 return(NULL); /* pro forma return */
611 }
612
613 /* Clear our BSDF representation and free memory */
614 void
615 clear_bsdf_rep(void)
616 {
617 while (mig_list != NULL) {
618 MIGRATION *mig = mig_list;
619 mig_list = mig->next;
620 free(mig);
621 }
622 while (dsf_list != NULL) {
623 RBFNODE *rbf = dsf_list;
624 dsf_list = rbf->next;
625 free(rbf);
626 }
627 bsdf_name[0] = '\0';
628 bsdf_manuf[0] = '\0';
629 inp_coverage = 0;
630 single_plane_incident = -1;
631 input_orient = output_orient = 0;
632 rbf_colorimetry = RBCunknown;
633 grid_res = GRIDRES;
634 memset(bsdf_hist, 0, sizeof(bsdf_hist));
635 bsdf_min = 0;
636 bsdf_spec_val = 0;
637 bsdf_spec_rad = 0;
638 }
639
640 /* Write our BSDF mesh interpolant out to the given binary stream */
641 void
642 save_bsdf_rep(FILE *ofp)
643 {
644 RBFNODE *rbf;
645 MIGRATION *mig;
646 int i, n;
647 /* finish header */
648 if (bsdf_name[0])
649 fprintf(ofp, "NAME=%s\n", bsdf_name);
650 if (bsdf_manuf[0])
651 fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf);
652 fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
653 fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
654 fprintf(ofp, "COLORIMETRY=%s\n", RBCident[rbf_colorimetry]);
655 fprintf(ofp, "GRIDRES=%d\n", grid_res);
656 fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
657 if ((bsdf_spec_val > bsdf_min) & (bsdf_spec_rad > 0))
658 fprintf(ofp, "BSDFSPEC= %f %f\n", bsdf_spec_val, bsdf_spec_rad);
659 fputformat(BSDFREP_FMT, ofp);
660 fputc('\n', ofp);
661 putint(BSDFREP_MAGIC, 2, ofp);
662 /* write each DSF */
663 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
664 putint(rbf->ord, 4, ofp);
665 putflt(rbf->invec[0], ofp);
666 putflt(rbf->invec[1], ofp);
667 putflt(rbf->invec[2], ofp);
668 putflt(rbf->vtotal, ofp);
669 putint(rbf->nrbf, 4, ofp);
670 for (i = 0; i < rbf->nrbf; i++) {
671 putflt(rbf->rbfa[i].peak, ofp);
672 putint(rbf->rbfa[i].chroma, 2, ofp);
673 putint(rbf->rbfa[i].crad, 2, ofp);
674 putint(rbf->rbfa[i].gx, 2, ofp);
675 putint(rbf->rbfa[i].gy, 2, ofp);
676 }
677 }
678 putint(-1, 4, ofp); /* terminator */
679 /* write each migration matrix */
680 for (mig = mig_list; mig != NULL; mig = mig->next) {
681 int zerocnt = 0;
682 putint(mig->rbfv[0]->ord, 4, ofp);
683 putint(mig->rbfv[1]->ord, 4, ofp);
684 /* write out as sparse data */
685 n = mtx_nrows(mig) * mtx_ncols(mig);
686 for (i = 0; i < n; i++) {
687 if (zerocnt == 0xff) {
688 putint(0xff, 1, ofp); zerocnt = 0;
689 }
690 if (mig->mtx[i] != 0) {
691 putint(zerocnt, 1, ofp); zerocnt = 0;
692 putflt(mig->mtx[i], ofp);
693 } else
694 ++zerocnt;
695 }
696 putint(zerocnt, 1, ofp);
697 }
698 putint(-1, 4, ofp); /* terminator */
699 putint(-1, 4, ofp);
700 if (fflush(ofp) == EOF) {
701 fprintf(stderr, "%s: error writing BSDF interpolant\n",
702 progname);
703 exit(1);
704 }
705 }
706
707 /* Check header line for critical information */
708 static int
709 headline(char *s, void *p)
710 {
711 char fmt[MAXFMTLEN];
712 int i;
713
714 if (isheadid(s))
715 return(0);
716 if (!strncmp(s, "NAME=", 5)) {
717 strcpy(bsdf_name, s+5);
718 bsdf_name[strlen(bsdf_name)-1] = '\0';
719 return(1);
720 }
721 if (!strncmp(s, "MANUFACT=", 9)) {
722 strcpy(bsdf_manuf, s+9);
723 bsdf_manuf[strlen(bsdf_manuf)-1] = '\0';
724 return(1);
725 }
726 if (!strncmp(s, "SYMMETRY=", 9)) {
727 inp_coverage = atoi(s+9);
728 single_plane_incident = !inp_coverage;
729 return(1);
730 }
731 if (!strncmp(s, "IO_SIDES=", 9)) {
732 sscanf(s+9, "%d %d", &input_orient, &output_orient);
733 return(1);
734 }
735 if (!strncmp(s, "COLORIMETRY=", 12)) {
736 fmt[0] = '\0';
737 sscanf(s+12, "%s", fmt);
738 for (i = RBCunknown; i >= 0; i--)
739 if (!strcmp(fmt, RBCident[i]))
740 break;
741 if (i < 0)
742 return(-1);
743 rbf_colorimetry = i;
744 return(1);
745 }
746 if (!strncmp(s, "GRIDRES=", 8)) {
747 sscanf(s+8, "%d", &grid_res);
748 return(1);
749 }
750 if (!strncmp(s, "BSDFMIN=", 8)) {
751 sscanf(s+8, "%lf", &bsdf_min);
752 return(1);
753 }
754 if (!strncmp(s, "BSDFSPEC=", 9)) {
755 sscanf(s+9, "%lf %lf", &bsdf_spec_val, &bsdf_spec_rad);
756 return(1);
757 }
758 if (formatval(fmt, s))
759 return (strcmp(fmt, BSDFREP_FMT) ? -1 : 0);
760 if (sir_headshare != NULL)
761 return ((*sir_headshare)(s));
762 return(0);
763 }
764
765 /* Read a BSDF mesh interpolant from the given binary stream */
766 int
767 load_bsdf_rep(FILE *ifp)
768 {
769 RBFNODE rbfh;
770 int from_ord, to_ord;
771 int i;
772
773 clear_bsdf_rep();
774 if (ifp == NULL)
775 return(0);
776 if (getheader(ifp, headline, NULL) < 0 || (single_plane_incident < 0) |
777 !input_orient | !output_orient |
778 (grid_res < 16) | (grid_res > 0xffff)) {
779 fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
780 progname);
781 return(0);
782 }
783 if (getint(2, ifp) != BSDFREP_MAGIC) {
784 fprintf(stderr, "%s: bad magic number for BSDF interpolant\n",
785 progname);
786 return(0);
787 }
788 memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */
789 while ((rbfh.ord = getint(4, ifp)) >= 0) {
790 RBFNODE *newrbf;
791
792 rbfh.invec[0] = getflt(ifp);
793 rbfh.invec[1] = getflt(ifp);
794 rbfh.invec[2] = getflt(ifp);
795 if (normalize(rbfh.invec) == 0) {
796 fprintf(stderr, "%s: zero incident vector\n", progname);
797 return(0);
798 }
799 rbfh.vtotal = getflt(ifp);
800 rbfh.nrbf = getint(4, ifp);
801 newrbf = (RBFNODE *)malloc(sizeof(RBFNODE) +
802 sizeof(RBFVAL)*(rbfh.nrbf-1));
803 if (newrbf == NULL)
804 goto memerr;
805 *newrbf = rbfh;
806 for (i = 0; i < rbfh.nrbf; i++) {
807 newrbf->rbfa[i].peak = getflt(ifp);
808 newrbf->rbfa[i].chroma = getint(2, ifp) & 0xffff;
809 newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
810 newrbf->rbfa[i].gx = getint(2, ifp) & 0xffff;
811 newrbf->rbfa[i].gy = getint(2, ifp) & 0xffff;
812 }
813 if (feof(ifp))
814 goto badEOF;
815 /* insert in global list */
816 if (insert_dsf(newrbf) != rbfh.ord) {
817 fprintf(stderr, "%s: error adding DSF\n", progname);
818 return(0);
819 }
820 }
821 /* read each migration matrix */
822 while ((from_ord = getint(4, ifp)) >= 0 &&
823 (to_ord = getint(4, ifp)) >= 0) {
824 RBFNODE *from_rbf = get_dsf(from_ord);
825 RBFNODE *to_rbf = get_dsf(to_ord);
826 MIGRATION *newmig;
827 int n;
828
829 if ((from_rbf == NULL) | (to_rbf == NULL)) {
830 fprintf(stderr,
831 "%s: bad DSF reference in migration edge\n",
832 progname);
833 return(0);
834 }
835 n = from_rbf->nrbf * to_rbf->nrbf;
836 newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
837 sizeof(float)*(n-1));
838 if (newmig == NULL)
839 goto memerr;
840 newmig->rbfv[0] = from_rbf;
841 newmig->rbfv[1] = to_rbf;
842 memset(newmig->mtx, 0, sizeof(float)*n);
843 for (i = 0; ; ) { /* read sparse data */
844 int zc = getint(1, ifp) & 0xff;
845 if ((i += zc) >= n)
846 break;
847 if (zc == 0xff)
848 continue;
849 newmig->mtx[i++] = getflt(ifp);
850 }
851 if (feof(ifp))
852 goto badEOF;
853 /* insert in edge lists */
854 newmig->enxt[0] = from_rbf->ejl;
855 from_rbf->ejl = newmig;
856 newmig->enxt[1] = to_rbf->ejl;
857 to_rbf->ejl = newmig;
858 /* push onto global list */
859 newmig->next = mig_list;
860 mig_list = newmig;
861 }
862 return(1); /* success! */
863 memerr:
864 fprintf(stderr, "%s: Out of memory in load_bsdf_rep()\n", progname);
865 exit(1);
866 badEOF:
867 fprintf(stderr, "%s: Unexpected EOF in load_bsdf_rep()\n", progname);
868 return(0);
869 }