ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfrep.c
(Generate patch)

Comparing ray/src/cv/bsdfrep.c (file contents):
Revision 2.8 by greg, Fri Nov 9 02:16:29 2012 UTC vs.
Revision 2.38 by greg, Sun May 18 01:46:05 2025 UTC

# Line 9 | Line 9 | static const char RCSid[] = "$Id$";
9  
10   #define _USE_MATH_DEFINES
11   #include <stdlib.h>
12 – #include <string.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  
# Line 26 | Line 29 | int                    single_plane_incident = -1;
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  
# Line 35 | Line 53 | MIGRATION              *mig_list = NULL;
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   {
42 –        if (!input_orient)              /* check input orientation */
43 –                input_orient = 1 - 2*(new_theta > 90.);
44 –        else if (input_orient > 0 ^ new_theta < 90.) {
45 –                fprintf(stderr,
46 –                "%s: Cannot handle input angles on both sides of surface\n",
47 –                                progname);
48 –                return(0);
49 –        }
63                                          /* normalize angle ranges */
64          while (new_theta < -180.)
65                  new_theta += 360.;
# Line 56 | Line 69 | new_input_direction(double new_theta, double new_phi)
69                  new_theta = -new_theta;
70                  new_phi += 180.;
71          }
59 –        if ((theta_in_deg = new_theta) < 1.0)
60 –                return(1);              /* don't rely on phi near normal */
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)
# Line 83 | Line 105 | 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:
# Line 189 | Line 214 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
214          int                     pos[2];
215          int                     n;
216  
217 <        for (n = rbf->nrbf; n-- > 0; ) {
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);
# Line 199 | Line 224 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
224          VCOPY(rbf->invec, invec);
225   }
226  
227 < /* Compute volume associated with Gaussian lobe */
228 < double
229 < rbf_volume(const RBFVAL *rbfp)
230 < {
231 <        double  rad = R2ANG(rbfp->crad);
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 <        return((2.*M_PI) * rbfp->peak * rad*rad);
236 < }
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 < /* Compute outgoing vector from grid position */
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 < {
271 <        double  uv[2];
270 > {                               /* table-free version */
271 >        RREAL   uv[2];
272          double  r2;
273 <        
274 <        SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
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);
# Line 223 | Line 282 | ovec_from_pos(FVECT vec, int xpos, int ypos)
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 <        double  sq[2];          /* uniform hemispherical projection */
291 >        RREAL   sq[2];          /* uniform hemispherical projection */
292          double  norm = 1./sqrt(1. + fabs(vec[2]));
293  
294 <        SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
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 < /* Evaluate RBF for DSF at the given normalized outgoing direction */
300 > /* Compute volume associated with Gaussian lobe */
301   double
302 < eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
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          sig2;
330 >        double          rad2;
331          int             n;
332 <
333 <        if (rp == NULL)
334 <                return(.0);
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 <                sig2 = R2ANG(rbfp->crad);
356 <                sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
357 <                if (sig2 > -19.)
358 <                        res += rbfp->peak * exp(sig2);
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 <        return(res);
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)
# Line 270 | Line 392 | insert_dsf(RBFNODE *newrbf)
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)\n",
396 <                                        progname);
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                  }
# Line 359 | Line 482 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
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)
# Line 373 | Line 624 | clear_bsdf_rep(void)
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 */
# Line 387 | Line 645 | save_bsdf_rep(FILE *ofp)
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);
# Line 402 | Line 669 | save_bsdf_rep(FILE *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, 1, ofp);
675 <                        putint(rbf->rbfa[i].gy, 1, 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 */
# Line 440 | Line 708 | save_bsdf_rep(FILE *ofp)
708   static int
709   headline(char *s, void *p)
710   {
711 <        char    fmt[32];
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(0);
729 >                return(1);
730          }
731          if (!strncmp(s, "IO_SIDES=", 9)) {
732                  sscanf(s+9, "%d %d", &input_orient, &output_orient);
733 <                return(0);
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(0);
748 >                return(1);
749          }
750 <        if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
751 <                return(-1);
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  
# Line 471 | Line 773 | load_bsdf_rep(FILE *ifp)
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) {
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 <        rbfh.next = NULL;               /* read each DSF */
784 <        rbfh.ejl = NULL;
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 <                memcpy(newrbf, &rbfh, sizeof(RBFNODE)-sizeof(RBFVAL));
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(1, ifp) & 0xff;
811 <                        newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
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;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines