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

Comparing ray/src/cv/bsdfrep.c (file contents):
Revision 2.21 by greg, Fri Mar 7 21:21:02 2014 UTC vs.
Revision 2.33 by greg, Sat Dec 28 18:05:14 2019 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 48 | Line 57 | double                 theta_in_deg, phi_in_deg;
57   int
58   new_input_direction(double new_theta, double new_phi)
59   {
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 –        }
60                                          /* normalize angle ranges */
61          while (new_theta < -180.)
62                  new_theta += 360.;
# Line 65 | Line 66 | new_input_direction(double new_theta, double new_phi)
66                  new_theta = -new_theta;
67                  new_phi += 180.;
68          }
68 –        if ((theta_in_deg = new_theta) < 1.0)
69 –                return(1);              /* don't rely on phi near normal */
69          while (new_phi < 0)
70                  new_phi += 360.;
71          while (new_phi >= 360.)
72                  new_phi -= 360.;
73 +                                        /* check input orientation */
74 +        if (!input_orient)
75 +                input_orient = 1 - 2*(new_theta > 90.);
76 +        else if (input_orient > 0 ^ new_theta < 90.) {
77 +                fprintf(stderr,
78 +                "%s: Cannot handle input angles on both sides of surface\n",
79 +                                progname);
80 +                return(0);
81 +        }
82 +        if ((theta_in_deg = new_theta) < 1.0)
83 +                return(1);              /* don't rely on phi near normal */
84          if (single_plane_incident > 0)  /* check input coverage */
85                  single_plane_incident = (round(new_phi) == round(phi_in_deg));
86          else if (single_plane_incident < 0)
# Line 198 | Line 208 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
208          int                     pos[2];
209          int                     n;
210  
211 <        for (n = ((-.01 > phi) | (phi > .01))*rbf->nrbf; n-- > 0; ) {
211 >        for (n = (cos(phi) < 1.-FTINY)*rbf->nrbf; n-- > 0; ) {
212                  ovec_from_pos(outvec, rbf->rbfa[n].gx, rbf->rbfa[n].gy);
213                  spinvector(outvec, outvec, vnorm, phi);
214                  pos_from_vec(pos, outvec);
# Line 257 | Line 267 | rbf_volume(const RBFVAL *rbfp)
267          return(integ);
268   }
269  
270 < /* Evaluate RBF for DSF at the given normalized outgoing direction */
271 < double
272 < eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
270 > /* Evaluate BSDF at the given normalized outgoing direction in color */
271 > SDError
272 > eval_rbfcol(SDValue *sv, const RBFNODE *rp, const FVECT outvec)
273   {
274          const double    rfact2 = (38./M_PI/M_PI)*(grid_res*grid_res);
265 –        double          minval = bsdf_min*output_orient*outvec[2];
275          int             pos[2];
276          double          res = 0;
277 +        double          usum = 0, vsum = 0;
278          const RBFVAL    *rbfp;
279          FVECT           odir;
280          double          rad2;
281          int             n;
282 +                                /* assign default value */
283 +        sv->spec = c_dfcolor;
284 +        sv->cieY = bsdf_min;
285                                  /* check for wrong side */
286 <        if (outvec[2] > 0 ^ output_orient > 0)
287 <                return(.0);
288 <                                /* use minimum if no information avail. */
289 <        if (rp == NULL)
290 <                return(minval);
286 >        if (outvec[2] > 0 ^ output_orient > 0) {
287 >                strcpy(SDerrorDetail, "Wrong-side scattering query");
288 >                return(SDEargument);
289 >        }
290 >        if (rp == NULL)         /* return minimum if no information avail. */
291 >                return(SDEnone);
292                                  /* optimization for fast lobe culling */
293          pos_from_vec(pos, outvec);
294                                  /* sum radial basis function */
# Line 282 | Line 296 | eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
296          for (n = rp->nrbf; n--; rbfp++) {
297                  int     d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
298                                  (pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
299 +                double  val;
300                  rad2 = R2ANG(rbfp->crad);
301                  rad2 *= rad2;
302                  if (d2 > rad2*rfact2)
303                          continue;
304                  ovec_from_pos(odir, rbfp->gx, rbfp->gy);
305 <                res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
305 >                val = rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
306 >                if (rbf_colorimetry == RBCtristimulus) {
307 >                        usum += val * (rbfp->chroma & 0xff);
308 >                        vsum += val * (rbfp->chroma>>8 & 0xff);
309 >                }
310 >                res += val;
311          }
312 <        if (res < minval)       /* never return less than minval */
313 <                return(minval);
314 <        return(res);
312 >        sv->cieY = res / COSF(outvec[2]);
313 >        if (sv->cieY < bsdf_min) {      /* never return less than bsdf_min */
314 >                sv->cieY = bsdf_min;
315 >        } else if (rbf_colorimetry == RBCtristimulus) {
316 >                C_CHROMA        cres = (int)(usum/res + frandom());
317 >                cres |= (int)(vsum/res + frandom()) << 8;
318 >                c_decodeChroma(&sv->spec, cres);
319 >        }
320 >        return(SDEnone);
321   }
322  
323 + /* Evaluate BSDF at the given normalized outgoing direction in Y */
324 + double
325 + eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
326 + {
327 +        SDValue sv;
328 +
329 +        if (eval_rbfcol(&sv, rp, outvec) == SDEnone)
330 +                return(sv.cieY);
331 +
332 +        return(0.0);
333 + }
334 +
335   /* Insert a new directional scattering function in our global list */
336   int
337   insert_dsf(RBFNODE *newrbf)
# Line 304 | Line 342 | insert_dsf(RBFNODE *newrbf)
342          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
343                  if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
344                          fprintf(stderr,
345 <                                "%s: Duplicate incident measurement (ignored)\n",
346 <                                        progname);
345 >                "%s: Duplicate incident measurement ignored at (%.1f,%.1f)\n",
346 >                                        progname, get_theta180(newrbf->invec),
347 >                                        get_phi360(newrbf->invec));
348                          free(newrbf);
349                          return(-1);
350                  }
# Line 393 | Line 432 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
432          return((rbfv[0] != NULL) + (rbfv[1] != NULL));
433   }
434  
435 + /* Return single-lobe specular RBF for the given incident direction */
436 + RBFNODE *
437 + def_rbf_spec(const FVECT invec)
438 + {
439 +        RBFNODE         *rbf;
440 +        FVECT           ovec;
441 +        int             pos[2];
442 +
443 +        if (input_orient > 0 ^ invec[2] > 0)    /* wrong side? */
444 +                return(NULL);
445 +        if ((bsdf_spec_val <= bsdf_min) | (bsdf_spec_rad <= 0))
446 +                return(NULL);                   /* nothing set */
447 +        rbf = (RBFNODE *)malloc(sizeof(RBFNODE));
448 +        if (rbf == NULL)
449 +                return(NULL);
450 +        ovec[0] = -invec[0];
451 +        ovec[1] = -invec[1];
452 +        ovec[2] = invec[2]*(2*(input_orient==output_orient) - 1);
453 +        pos_from_vec(pos, ovec);
454 +        rbf->ord = 0;
455 +        rbf->next = NULL;
456 +        rbf->ejl = NULL;
457 +        VCOPY(rbf->invec, invec);
458 +        rbf->nrbf = 1;
459 +        rbf->rbfa[0].peak = bsdf_spec_val * COSF(ovec[2]);
460 +        rbf->rbfa[0].chroma = c_dfchroma;
461 +        rbf->rbfa[0].crad = ANG2R(bsdf_spec_rad);
462 +        rbf->rbfa[0].gx = pos[0];
463 +        rbf->rbfa[0].gy = pos[1];
464 +        rbf->vtotal = rbf_volume(rbf->rbfa);
465 +        return(rbf);
466 + }
467 +
468   /* Advect and allocate new RBF along edge (internal call) */
469   RBFNODE *
470   e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
# Line 449 | Line 521 | tryagain:
521              const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
522              const float         peak0 = rbf0i->peak;
523              const double        rad0 = R2ANG(rbf0i->crad);
524 +            C_COLOR             cc0;
525              FVECT               v0;
526              float               mv;
527              ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
528 +            c_decodeChroma(&cc0, rbf0i->chroma);
529              for (j = 0; j < mtx_ncols(mig); j++)
530                  if ((mv = mtx_coef(mig,i,j)) > cthresh) {
531                          const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
# Line 462 | Line 536 | tryagain:
536                          rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
537                          rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
538                                                  rad0*rad0/rad2;
539 +                        if (rbf_colorimetry == RBCtristimulus) {
540 +                                C_COLOR cres;
541 +                                c_decodeChroma(&cres, rbf1j->chroma);
542 +                                c_cmix(&cres, 1.-t, &cc0, t, &cres);
543 +                                rbf->rbfa[n].chroma = c_encodeChroma(&cres);
544 +                        } else
545 +                                rbf->rbfa[n].chroma = c_dfchroma;
546                          rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
547                          ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
548                          geodesic(v, v0, v, t, GEOD_REL);
# Line 498 | Line 579 | clear_bsdf_rep(void)
579          inp_coverage = 0;
580          single_plane_incident = -1;
581          input_orient = output_orient = 0;
582 +        rbf_colorimetry = RBCunknown;
583          grid_res = GRIDRES;
584 +        memset(bsdf_hist, 0, sizeof(bsdf_hist));
585 +        bsdf_min = 0;
586 +        bsdf_spec_val = 0;
587 +        bsdf_spec_rad = 0;
588   }
589  
590   /* Write our BSDF mesh interpolant out to the given binary stream */
# Line 515 | Line 601 | save_bsdf_rep(FILE *ofp)
601                  fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf);
602          fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
603          fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
604 +        fprintf(ofp, "COLORIMETRY=%s\n", RBCident[rbf_colorimetry]);
605          fprintf(ofp, "GRIDRES=%d\n", grid_res);
606          fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
607 +        if ((bsdf_spec_val > bsdf_min) & (bsdf_spec_rad > 0))
608 +                fprintf(ofp, "BSDFSPEC= %f %f\n", bsdf_spec_val, bsdf_spec_rad);
609          fputformat(BSDFREP_FMT, ofp);
610          fputc('\n', ofp);
611 +        putint(BSDFREP_MAGIC, 2, ofp);
612                                          /* write each DSF */
613          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
614                  putint(rbf->ord, 4, ofp);
# Line 529 | Line 619 | save_bsdf_rep(FILE *ofp)
619                  putint(rbf->nrbf, 4, ofp);
620                  for (i = 0; i < rbf->nrbf; i++) {
621                          putflt(rbf->rbfa[i].peak, ofp);
622 +                        putint(rbf->rbfa[i].chroma, 2, ofp);
623                          putint(rbf->rbfa[i].crad, 2, ofp);
624 <                        putint(rbf->rbfa[i].gx, 1, ofp);
625 <                        putint(rbf->rbfa[i].gy, 1, ofp);
624 >                        putint(rbf->rbfa[i].gx, 2, ofp);
625 >                        putint(rbf->rbfa[i].gy, 2, ofp);
626                  }
627          }
628          putint(-1, 4, ofp);             /* terminator */
# Line 567 | Line 658 | save_bsdf_rep(FILE *ofp)
658   static int
659   headline(char *s, void *p)
660   {
661 <        char    fmt[32];
661 >        char    fmt[MAXFMTLEN];
662 >        int     i;
663  
664          if (!strncmp(s, "NAME=", 5)) {
665                  strcpy(bsdf_name, s+5);
# Line 586 | Line 678 | headline(char *s, void *p)
678                  sscanf(s+9, "%d %d", &input_orient, &output_orient);
679                  return(0);
680          }
681 +        if (!strncmp(s, "COLORIMETRY=", 12)) {
682 +                fmt[0] = '\0';
683 +                sscanf(s+12, "%s", fmt);
684 +                for (i = RBCunknown; i >= 0; i--)
685 +                        if (!strcmp(fmt, RBCident[i]))
686 +                                break;
687 +                if (i < 0)
688 +                        return(-1);
689 +                rbf_colorimetry = i;
690 +                return(0);
691 +        }
692          if (!strncmp(s, "GRIDRES=", 8)) {
693                  sscanf(s+8, "%d", &grid_res);
694                  return(0);
# Line 594 | Line 697 | headline(char *s, void *p)
697                  sscanf(s+8, "%lf", &bsdf_min);
698                  return(0);
699          }
700 +        if (!strncmp(s, "BSDFSPEC=", 9)) {
701 +                sscanf(s+9, "%lf %lf", &bsdf_spec_val, &bsdf_spec_rad);
702 +                return(0);
703 +        }
704          if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
705                  return(-1);
706          return(0);
# Line 611 | Line 718 | load_bsdf_rep(FILE *ifp)
718          if (ifp == NULL)
719                  return(0);
720          if (getheader(ifp, headline, NULL) < 0 || (single_plane_incident < 0) |
721 <                        !input_orient | !output_orient) {
721 >                        !input_orient | !output_orient |
722 >                        (grid_res < 16) | (grid_res > 0xffff)) {
723                  fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
724                                  progname);
725                  return(0);
726          }
727 +        if (getint(2, ifp) != BSDFREP_MAGIC) {
728 +                fprintf(stderr, "%s: bad magic number for BSDF interpolant\n",
729 +                                progname);
730 +                return(0);
731 +        }
732          memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */
733          while ((rbfh.ord = getint(4, ifp)) >= 0) {
734                  RBFNODE         *newrbf;
# Line 636 | Line 749 | load_bsdf_rep(FILE *ifp)
749                  *newrbf = rbfh;
750                  for (i = 0; i < rbfh.nrbf; i++) {
751                          newrbf->rbfa[i].peak = getflt(ifp);
752 +                        newrbf->rbfa[i].chroma = getint(2, ifp) & 0xffff;
753                          newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
754 <                        newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
755 <                        newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
754 >                        newrbf->rbfa[i].gx = getint(2, ifp) & 0xffff;
755 >                        newrbf->rbfa[i].gy = getint(2, ifp) & 0xffff;
756                  }
757                  if (feof(ifp))
758                          goto badEOF;

Diff Legend

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