ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfrep.c
Revision: 2.31
Committed: Wed Feb 3 18:53:14 2016 UTC (8 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R1
Changes since 2.30: +6 -5 lines
Log Message:
Return minimum BSDF value as gray

File Contents

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