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.17 by greg, Thu Oct 24 16:11:38 2013 UTC vs.
Revision 2.31 by greg, Wed Feb 3 18:53:14 2016 UTC

# Line 14 | Line 14 | static const char RCSid[] = "$Id$";
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  
# Line 26 | Line 30 | int                    single_plane_incident = -1;
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;
# Line 45 | Line 58 | double                 theta_in_deg, phi_in_deg;
58   int
59   new_input_direction(double new_theta, double new_phi)
60   {
48 –        if (!input_orient)              /* check input orientation */
49 –                input_orient = 1 - 2*(new_theta > 90.);
50 –        else if (input_orient > 0 ^ new_theta < 90.) {
51 –                fprintf(stderr,
52 –                "%s: Cannot handle input angles on both sides of surface\n",
53 –                                progname);
54 –                return(0);
55 –        }
61                                          /* normalize angle ranges */
62          while (new_theta < -180.)
63                  new_theta += 360.;
# Line 62 | Line 67 | new_input_direction(double new_theta, double new_phi)
67                  new_theta = -new_theta;
68                  new_phi += 180.;
69          }
65 –        if ((theta_in_deg = new_theta) < 1.0)
66 –                return(1);              /* don't rely on phi near normal */
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)
# Line 195 | Line 209 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
209          int                     pos[2];
210          int                     n;
211  
212 <        for (n = ((-.01 > phi) | (phi > .01))*rbf->nrbf; n-- > 0; ) {
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);
# Line 254 | Line 268 | rbf_volume(const RBFVAL *rbfp)
268          return(integ);
269   }
270  
271 < /* Evaluate RBF for DSF at the given normalized outgoing direction */
272 < double
273 < eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
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);
262 –        double          minval = bsdf_min*output_orient*outvec[2];
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 <                return(.0);
289 <                                /* use minimum if no information avail. */
290 <        if (rp == NULL)
291 <                return(minval);
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 */
# Line 279 | Line 297 | eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
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 <                res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
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 <        if (res < minval)       /* never return less than minval */
314 <                return(minval);
315 <        return(res);
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)
# Line 301 | Line 343 | insert_dsf(RBFNODE *newrbf)
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)\n",
347 <                                        progname);
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                  }
# Line 390 | Line 433 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
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)
# Line 404 | Line 575 | clear_bsdf_rep(void)
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 */
# Line 418 | Line 596 | save_bsdf_rep(FILE *ofp)
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);
# Line 434 | Line 620 | save_bsdf_rep(FILE *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, 1, ofp);
626 <                        putint(rbf->rbfa[i].gy, 1, 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 */
# Line 472 | Line 659 | save_bsdf_rep(FILE *ofp)
659   static int
660   headline(char *s, void *p)
661   {
662 <        char    fmt[32];
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;
# Line 483 | Line 679 | headline(char *s, void *p)
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);
# Line 491 | Line 698 | headline(char *s, void *p)
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);
# Line 507 | Line 718 | load_bsdf_rep(FILE *ifp)
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) {
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 <        rbfh.next = NULL;               /* read each DSF */
729 <        rbfh.ejl = NULL;
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  
# Line 531 | Line 747 | load_bsdf_rep(FILE *ifp)
747                                          sizeof(RBFVAL)*(rbfh.nrbf-1));
748                  if (newrbf == NULL)
749                          goto memerr;
750 <                memcpy(newrbf, &rbfh, sizeof(RBFNODE)-sizeof(RBFVAL));
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(1, ifp) & 0xff;
756 <                        newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
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;

Diff Legend

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