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.24 by greg, Tue Mar 25 14:55:35 2014 UTC vs.
Revision 2.37 by greg, Wed Dec 15 01:38:50 2021 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];
# Line 29 | 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;
# Line 44 | 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   {
51 –        if (!input_orient)              /* check input orientation */
52 –                input_orient = 1 - 2*(new_theta > 90.);
53 –        else if (input_orient > 0 ^ new_theta < 90.) {
54 –                fprintf(stderr,
55 –                "%s: Cannot handle input angles on both sides of surface\n",
56 –                                progname);
57 –                return(0);
58 –        }
63                                          /* normalize angle ranges */
64          while (new_theta < -180.)
65                  new_theta += 360.;
# Line 65 | Line 69 | new_input_direction(double new_theta, double new_phi)
69                  new_theta = -new_theta;
70                  new_phi += 180.;
71          }
68 –        if ((theta_in_deg = new_theta) < 1.0)
69 –                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 91 | Line 104 | new_input_direction(double new_theta, double new_phi)
104   int
105   use_symmetry(FVECT vec)
106   {
107 <        const double    phi = get_phi360(vec);
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 212 | Line 228 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
228   void
229   ovec_from_pos(FVECT vec, int xpos, int ypos)
230   {
231 <        double  uv[2];
231 >        RREAL   uv[2];
232          double  r2;
233          
234 <        SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
234 >        square2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
235                                  /* uniform hemispherical projection */
236          r2 = uv[0]*uv[0] + uv[1]*uv[1];
237          vec[0] = vec[1] = sqrt(2. - r2);
# Line 228 | Line 244 | ovec_from_pos(FVECT vec, int xpos, int ypos)
244   void
245   pos_from_vec(int pos[2], const FVECT vec)
246   {
247 <        double  sq[2];          /* uniform hemispherical projection */
247 >        RREAL   sq[2];          /* uniform hemispherical projection */
248          double  norm = 1./sqrt(1. + fabs(vec[2]));
249  
250 <        SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
250 >        disk2square(sq, vec[0]*norm, vec[1]*norm);
251  
252          pos[0] = (int)(sq[0]*grid_res);
253          pos[1] = (int)(sq[1]*grid_res);
# Line 257 | Line 273 | rbf_volume(const RBFVAL *rbfp)
273          return(integ);
274   }
275  
276 < /* Evaluate RBF for DSF at the given normalized outgoing direction */
277 < double
278 < eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
276 > /* Evaluate BSDF at the given normalized outgoing direction in color */
277 > SDError
278 > eval_rbfcol(SDValue *sv, const RBFNODE *rp, const FVECT outvec)
279   {
280          const double    rfact2 = (38./M_PI/M_PI)*(grid_res*grid_res);
265 –        double          minval = bsdf_min*output_orient*outvec[2];
281          int             pos[2];
282          double          res = 0;
283 +        double          usum = 0, vsum = 0;
284          const RBFVAL    *rbfp;
285          FVECT           odir;
286          double          rad2;
287          int             n;
288 +                                /* assign default value */
289 +        sv->spec = c_dfcolor;
290 +        sv->cieY = bsdf_min;
291                                  /* check for wrong side */
292 <        if (outvec[2] > 0 ^ output_orient > 0)
293 <                return(.0);
294 <                                /* use minimum if no information avail. */
295 <        if (rp == NULL)
296 <                return(minval);
292 >        if (outvec[2] > 0 ^ output_orient > 0) {
293 >                strcpy(SDerrorDetail, "Wrong-side scattering query");
294 >                return(SDEargument);
295 >        }
296 >        if (rp == NULL)         /* return minimum if no information avail. */
297 >                return(SDEnone);
298                                  /* optimization for fast lobe culling */
299          pos_from_vec(pos, outvec);
300                                  /* sum radial basis function */
# Line 282 | Line 302 | eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
302          for (n = rp->nrbf; n--; rbfp++) {
303                  int     d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
304                                  (pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
305 +                double  val;
306                  rad2 = R2ANG(rbfp->crad);
307                  rad2 *= rad2;
308                  if (d2 > rad2*rfact2)
309                          continue;
310                  ovec_from_pos(odir, rbfp->gx, rbfp->gy);
311 <                res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
311 >                val = rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
312 >                if (rbf_colorimetry == RBCtristimulus) {
313 >                        usum += val * (rbfp->chroma & 0xff);
314 >                        vsum += val * (rbfp->chroma>>8 & 0xff);
315 >                }
316 >                res += val;
317          }
318 <        if (res < minval)       /* never return less than minval */
319 <                return(minval);
320 <        return(res);
318 >        sv->cieY = res / COSF(outvec[2]);
319 >        if (sv->cieY < bsdf_min) {      /* never return less than bsdf_min */
320 >                sv->cieY = bsdf_min;
321 >        } else if (rbf_colorimetry == RBCtristimulus) {
322 >                C_CHROMA        cres = (int)(usum/res + frandom());
323 >                cres |= (int)(vsum/res + frandom()) << 8;
324 >                c_decodeChroma(&sv->spec, cres);
325 >        }
326 >        return(SDEnone);
327   }
328  
329 + /* Evaluate BSDF at the given normalized outgoing direction in Y */
330 + double
331 + eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
332 + {
333 +        SDValue sv;
334 +
335 +        if (eval_rbfcol(&sv, rp, outvec) == SDEnone)
336 +                return(sv.cieY);
337 +
338 +        return(0.0);
339 + }
340 +
341   /* Insert a new directional scattering function in our global list */
342   int
343   insert_dsf(RBFNODE *newrbf)
# Line 394 | Line 438 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
438          return((rbfv[0] != NULL) + (rbfv[1] != NULL));
439   }
440  
441 + /* Return single-lobe specular RBF for the given incident direction */
442 + RBFNODE *
443 + def_rbf_spec(const FVECT invec)
444 + {
445 +        RBFNODE         *rbf;
446 +        FVECT           ovec;
447 +        int             pos[2];
448 +
449 +        if (input_orient > 0 ^ invec[2] > 0)    /* wrong side? */
450 +                return(NULL);
451 +        if ((bsdf_spec_val <= bsdf_min) | (bsdf_spec_rad <= 0))
452 +                return(NULL);                   /* nothing set */
453 +        rbf = (RBFNODE *)malloc(sizeof(RBFNODE));
454 +        if (rbf == NULL)
455 +                return(NULL);
456 +        ovec[0] = -invec[0];
457 +        ovec[1] = -invec[1];
458 +        ovec[2] = invec[2]*(2*(input_orient==output_orient) - 1);
459 +        pos_from_vec(pos, ovec);
460 +        rbf->ord = 0;
461 +        rbf->next = NULL;
462 +        rbf->ejl = NULL;
463 +        VCOPY(rbf->invec, invec);
464 +        rbf->nrbf = 1;
465 +        rbf->rbfa[0].peak = bsdf_spec_val * COSF(ovec[2]);
466 +        rbf->rbfa[0].chroma = c_dfchroma;
467 +        rbf->rbfa[0].crad = ANG2R(bsdf_spec_rad);
468 +        rbf->rbfa[0].gx = pos[0];
469 +        rbf->rbfa[0].gy = pos[1];
470 +        rbf->vtotal = rbf_volume(rbf->rbfa);
471 +        return(rbf);
472 + }
473 +
474   /* Advect and allocate new RBF along edge (internal call) */
475   RBFNODE *
476   e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
# Line 404 | Line 481 | e_advect_rbf(const MIGRATION *mig, const FVECT invec,
481          double          t, full_dist;
482                                                  /* get relative position */
483          t = Acos(DOT(invec, mig->rbfv[0]->invec));
484 <        if (t < M_PI/grid_res) {                /* near first DSF */
484 >        if (t <= .001) {                        /* near first DSF */
485                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
486                  rbf = (RBFNODE *)malloc(n);
487                  if (rbf == NULL)
# Line 414 | Line 491 | e_advect_rbf(const MIGRATION *mig, const FVECT invec,
491                  return(rbf);
492          }
493          full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
494 <        if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
494 >        if (t >= full_dist-.001) {              /* near second DSF */
495                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
496                  rbf = (RBFNODE *)malloc(n);
497                  if (rbf == NULL)
# Line 450 | Line 527 | tryagain:
527              const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
528              const float         peak0 = rbf0i->peak;
529              const double        rad0 = R2ANG(rbf0i->crad);
530 +            C_COLOR             cc0;
531              FVECT               v0;
532              float               mv;
533              ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
534 +            c_decodeChroma(&cc0, rbf0i->chroma);
535              for (j = 0; j < mtx_ncols(mig); j++)
536                  if ((mv = mtx_coef(mig,i,j)) > cthresh) {
537                          const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
# Line 463 | Line 542 | tryagain:
542                          rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
543                          rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
544                                                  rad0*rad0/rad2;
545 +                        if (rbf_colorimetry == RBCtristimulus) {
546 +                                C_COLOR cres;
547 +                                c_decodeChroma(&cres, rbf1j->chroma);
548 +                                c_cmix(&cres, 1.-t, &cc0, t, &cres);
549 +                                rbf->rbfa[n].chroma = c_encodeChroma(&cres);
550 +                        } else
551 +                                rbf->rbfa[n].chroma = c_dfchroma;
552                          rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
553                          ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
554                          geodesic(v, v0, v, t, GEOD_REL);
# Line 499 | Line 585 | clear_bsdf_rep(void)
585          inp_coverage = 0;
586          single_plane_incident = -1;
587          input_orient = output_orient = 0;
588 +        rbf_colorimetry = RBCunknown;
589          grid_res = GRIDRES;
590 +        memset(bsdf_hist, 0, sizeof(bsdf_hist));
591 +        bsdf_min = 0;
592 +        bsdf_spec_val = 0;
593 +        bsdf_spec_rad = 0;
594   }
595  
596   /* Write our BSDF mesh interpolant out to the given binary stream */
# Line 516 | Line 607 | save_bsdf_rep(FILE *ofp)
607                  fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf);
608          fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
609          fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
610 +        fprintf(ofp, "COLORIMETRY=%s\n", RBCident[rbf_colorimetry]);
611          fprintf(ofp, "GRIDRES=%d\n", grid_res);
612          fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
613 +        if ((bsdf_spec_val > bsdf_min) & (bsdf_spec_rad > 0))
614 +                fprintf(ofp, "BSDFSPEC= %f %f\n", bsdf_spec_val, bsdf_spec_rad);
615          fputformat(BSDFREP_FMT, ofp);
616          fputc('\n', ofp);
617 +        putint(BSDFREP_MAGIC, 2, ofp);
618                                          /* write each DSF */
619          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
620                  putint(rbf->ord, 4, ofp);
# Line 530 | Line 625 | save_bsdf_rep(FILE *ofp)
625                  putint(rbf->nrbf, 4, ofp);
626                  for (i = 0; i < rbf->nrbf; i++) {
627                          putflt(rbf->rbfa[i].peak, ofp);
628 +                        putint(rbf->rbfa[i].chroma, 2, ofp);
629                          putint(rbf->rbfa[i].crad, 2, ofp);
630 <                        putint(rbf->rbfa[i].gx, 1, ofp);
631 <                        putint(rbf->rbfa[i].gy, 1, ofp);
630 >                        putint(rbf->rbfa[i].gx, 2, ofp);
631 >                        putint(rbf->rbfa[i].gy, 2, ofp);
632                  }
633          }
634          putint(-1, 4, ofp);             /* terminator */
# Line 568 | Line 664 | save_bsdf_rep(FILE *ofp)
664   static int
665   headline(char *s, void *p)
666   {
667 <        char    fmt[32];
667 >        char    fmt[MAXFMTLEN];
668 >        int     i;
669  
670 +        if (isheadid(s))
671 +                return(0);
672          if (!strncmp(s, "NAME=", 5)) {
673                  strcpy(bsdf_name, s+5);
674                  bsdf_name[strlen(bsdf_name)-1] = '\0';
675 +                return(1);
676          }
677          if (!strncmp(s, "MANUFACT=", 9)) {
678                  strcpy(bsdf_manuf, s+9);
679                  bsdf_manuf[strlen(bsdf_manuf)-1] = '\0';
680 +                return(1);
681          }
682          if (!strncmp(s, "SYMMETRY=", 9)) {
683                  inp_coverage = atoi(s+9);
684                  single_plane_incident = !inp_coverage;
685 <                return(0);
685 >                return(1);
686          }
687          if (!strncmp(s, "IO_SIDES=", 9)) {
688                  sscanf(s+9, "%d %d", &input_orient, &output_orient);
689 <                return(0);
689 >                return(1);
690          }
691 +        if (!strncmp(s, "COLORIMETRY=", 12)) {
692 +                fmt[0] = '\0';
693 +                sscanf(s+12, "%s", fmt);
694 +                for (i = RBCunknown; i >= 0; i--)
695 +                        if (!strcmp(fmt, RBCident[i]))
696 +                                break;
697 +                if (i < 0)
698 +                        return(-1);
699 +                rbf_colorimetry = i;
700 +                return(1);
701 +        }
702          if (!strncmp(s, "GRIDRES=", 8)) {
703                  sscanf(s+8, "%d", &grid_res);
704 <                return(0);
704 >                return(1);
705          }
706          if (!strncmp(s, "BSDFMIN=", 8)) {
707                  sscanf(s+8, "%lf", &bsdf_min);
708 <                return(0);
708 >                return(1);
709          }
710 <        if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
711 <                return(-1);
710 >        if (!strncmp(s, "BSDFSPEC=", 9)) {
711 >                sscanf(s+9, "%lf %lf", &bsdf_spec_val, &bsdf_spec_rad);
712 >                return(1);
713 >        }
714 >        if (formatval(fmt, s))
715 >                return (strcmp(fmt, BSDFREP_FMT) ? -1 : 0);
716 >        if (sir_headshare != NULL)
717 >                return ((*sir_headshare)(s));
718          return(0);
719   }
720  
# Line 613 | Line 731 | load_bsdf_rep(FILE *ifp)
731                  return(0);
732          if (getheader(ifp, headline, NULL) < 0 || (single_plane_incident < 0) |
733                          !input_orient | !output_orient |
734 <                        (grid_res < 16) | (grid_res > 256)) {
734 >                        (grid_res < 16) | (grid_res > 0xffff)) {
735                  fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
736                                  progname);
737                  return(0);
738          }
739 +        if (getint(2, ifp) != BSDFREP_MAGIC) {
740 +                fprintf(stderr, "%s: bad magic number for BSDF interpolant\n",
741 +                                progname);
742 +                return(0);
743 +        }
744          memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */
745          while ((rbfh.ord = getint(4, ifp)) >= 0) {
746                  RBFNODE         *newrbf;
# Line 638 | Line 761 | load_bsdf_rep(FILE *ifp)
761                  *newrbf = rbfh;
762                  for (i = 0; i < rbfh.nrbf; i++) {
763                          newrbf->rbfa[i].peak = getflt(ifp);
764 +                        newrbf->rbfa[i].chroma = getint(2, ifp) & 0xffff;
765                          newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
766 <                        newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
767 <                        newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
766 >                        newrbf->rbfa[i].gx = getint(2, ifp) & 0xffff;
767 >                        newrbf->rbfa[i].gy = getint(2, ifp) & 0xffff;
768                  }
769                  if (feof(ifp))
770                          goto badEOF;

Diff Legend

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