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.7 by greg, Thu Nov 8 22:05:04 2012 UTC vs.
Revision 2.26 by greg, Thu Aug 21 13:44:05 2014 UTC

# Line 14 | Line 14 | static const char RCSid[] = "$Id$";
14   #include "rtio.h"
15   #include "resolu.h"
16   #include "bsdfrep.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 +                                /* BSDF histogram */
33 + unsigned long           bsdf_hist[HISTLEN];
34 +
35 +                                /* BSDF value for boundary regions */
36 + double                  bsdf_min = 0;
37 + double                  bsdf_spec_peak = 0;
38 + double                  bsdf_spec_rad = 0;
39 +
40                                  /* processed incident DSF measurements */
41   RBFNODE                 *dsf_list = NULL;
42  
# Line 82 | Line 93 | new_input_direction(double new_theta, double new_phi)
93   int
94   use_symmetry(FVECT vec)
95   {
96 <        double  phi = get_phi360(vec);
96 >        const double    phi = get_phi360(vec);
97  
98          switch (inp_coverage) {
99          case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
# Line 188 | Line 199 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
199          FVECT                   outvec;
200          int                     pos[2];
201          int                     n;
202 < #ifdef DEBUG
203 <        double  tdiff = 180./M_PI*fabs(acos(invec[2])-acos(rbf->invec[2]));
193 <        if (tdiff >= 1.5)
194 <                fprintf(stderr,
195 <                        "%s: Warning - rotated theta differs by %.1f degrees\n",
196 <                                progname, tdiff);
197 < #endif
198 <        for (n = rbf->nrbf; n-- > 0; ) {
202 >
203 >        for (n = (cos(phi) < 1.-FTINY)*rbf->nrbf; n-- > 0; ) {
204                  ovec_from_pos(outvec, rbf->rbfa[n].gx, rbf->rbfa[n].gy);
205                  spinvector(outvec, outvec, vnorm, phi);
206                  pos_from_vec(pos, outvec);
# Line 205 | Line 210 | rotate_rbf(RBFNODE *rbf, const FVECT invec)
210          VCOPY(rbf->invec, invec);
211   }
212  
208 – /* Compute volume associated with Gaussian lobe */
209 – double
210 – rbf_volume(const RBFVAL *rbfp)
211 – {
212 –        double  rad = R2ANG(rbfp->crad);
213 –
214 –        return((2.*M_PI) * rbfp->peak * rad*rad);
215 – }
216 –
213   /* Compute outgoing vector from grid position */
214   void
215   ovec_from_pos(FVECT vec, int xpos, int ypos)
# Line 221 | Line 217 | ovec_from_pos(FVECT vec, int xpos, int ypos)
217          double  uv[2];
218          double  r2;
219          
220 <        SDsquare2disk(uv, (1./grid_res)*(xpos+.5), (1./grid_res)*(ypos+.5));
220 >        SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res);
221                                  /* uniform hemispherical projection */
222          r2 = uv[0]*uv[0] + uv[1]*uv[1];
223          vec[0] = vec[1] = sqrt(2. - r2);
# Line 243 | Line 239 | pos_from_vec(int pos[2], const FVECT vec)
239          pos[1] = (int)(sq[1]*grid_res);
240   }
241  
242 < /* Evaluate RBF for DSF at the given normalized outgoing direction */
242 > /* Compute volume associated with Gaussian lobe */
243   double
244 + rbf_volume(const RBFVAL *rbfp)
245 + {
246 +        double  rad = R2ANG(rbfp->crad);
247 +        FVECT   odir;
248 +        double  elev, integ;
249 +                                /* infinite integral approximation */
250 +        integ = (2.*M_PI) * rbfp->peak * rad*rad;
251 +                                /* check if we're near horizon */
252 +        ovec_from_pos(odir, rbfp->gx, rbfp->gy);
253 +        elev = output_orient*odir[2];
254 +                                /* apply cut-off correction if > 1% */
255 +        if (elev < 2.8*rad) {
256 +                /* elev = asin(elev);   /* this is so crude, anyway... */
257 +                integ *= 1. - .5*exp(-.5*elev*elev/(rad*rad));
258 +        }
259 +        return(integ);
260 + }
261 +
262 + /* Evaluate BSDF at the given normalized outgoing direction */
263 + double
264   eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
265   {
266 <        double          res = .0;
266 >        const double    rfact2 = (38./M_PI/M_PI)*(grid_res*grid_res);
267 >        int             pos[2];
268 >        double          res = 0;
269          const RBFVAL    *rbfp;
270          FVECT           odir;
271 <        double          sig2;
271 >        double          rad2;
272          int             n;
273 <
274 <        if (rp == NULL)
273 >                                /* check for wrong side */
274 >        if (outvec[2] > 0 ^ output_orient > 0)
275                  return(.0);
276 +                                /* use minimum if no information avail. */
277 +        if (rp == NULL)
278 +                return(bsdf_min);
279 +                                /* optimization for fast lobe culling */
280 +        pos_from_vec(pos, outvec);
281 +                                /* sum radial basis function */
282          rbfp = rp->rbfa;
283          for (n = rp->nrbf; n--; rbfp++) {
284 +                int     d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
285 +                                (pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
286 +                rad2 = R2ANG(rbfp->crad);
287 +                rad2 *= rad2;
288 +                if (d2 > rad2*rfact2)
289 +                        continue;
290                  ovec_from_pos(odir, rbfp->gx, rbfp->gy);
291 <                sig2 = R2ANG(rbfp->crad);
262 <                sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
263 <                if (sig2 > -19.)
264 <                        res += rbfp->peak * exp(sig2);
291 >                res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
292          }
293 +        res /= output_orient*outvec[2];
294 +        if (res < bsdf_min)     /* never return less than bsdf_min */
295 +                return(bsdf_min);
296          return(res);
297   }
298  
# Line 276 | Line 306 | insert_dsf(RBFNODE *newrbf)
306          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
307                  if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
308                          fprintf(stderr,
309 <                                "%s: Duplicate incident measurement (ignored)\n",
310 <                                        progname);
309 >                "%s: Duplicate incident measurement ignored at (%.1f,%.1f)\n",
310 >                                        progname, get_theta180(newrbf->invec),
311 >                                        get_phi360(newrbf->invec));
312                          free(newrbf);
313                          return(-1);
314                  }
# Line 365 | Line 396 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
396          return((rbfv[0] != NULL) + (rbfv[1] != NULL));
397   }
398  
399 + /* Return single-lobe specular RBF for the given incident direction */
400 + RBFNODE *
401 + def_rbf_spec(const FVECT invec)
402 + {
403 +        RBFNODE         *rbf;
404 +        FVECT           ovec;
405 +        int             pos[2];
406 +
407 +        if (input_orient > 0 ^ invec[2] > 0)    /* wrong side? */
408 +                return(NULL);
409 +        if ((bsdf_spec_peak <= bsdf_min) | (bsdf_spec_rad <= 0))
410 +                return(NULL);                   /* nothing set */
411 +        rbf = (RBFNODE *)malloc(sizeof(RBFNODE));
412 +        if (rbf == NULL)
413 +                return(NULL);
414 +        ovec[0] = -invec[0];
415 +        ovec[1] = -invec[1];
416 +        ovec[2] = invec[2]*(2*(input_orient==output_orient) - 1);
417 +        pos_from_vec(pos, ovec);
418 +        rbf->ord = 0;
419 +        rbf->next = NULL;
420 +        rbf->ejl = NULL;
421 +        VCOPY(rbf->invec, invec);
422 +        rbf->nrbf = 1;
423 +        rbf->rbfa[0].peak = bsdf_spec_peak * output_orient*ovec[2];
424 +        rbf->rbfa[0].crad = ANG2R(bsdf_spec_rad);
425 +        rbf->rbfa[0].gx = pos[0];
426 +        rbf->rbfa[0].gy = pos[1];
427 +        rbf->vtotal = rbf_volume(rbf->rbfa);
428 +        return(rbf);
429 + }
430 +
431 + /* Advect and allocate new RBF along edge (internal call) */
432 + RBFNODE *
433 + e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
434 + {
435 +        double          cthresh = FTINY;
436 +        RBFNODE         *rbf;
437 +        int             n, i, j;
438 +        double          t, full_dist;
439 +                                                /* get relative position */
440 +        t = Acos(DOT(invec, mig->rbfv[0]->invec));
441 +        if (t < M_PI/grid_res) {                /* near first DSF */
442 +                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
443 +                rbf = (RBFNODE *)malloc(n);
444 +                if (rbf == NULL)
445 +                        goto memerr;
446 +                memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
447 +                rbf->next = NULL; rbf->ejl = NULL;
448 +                return(rbf);
449 +        }
450 +        full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
451 +        if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
452 +                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
453 +                rbf = (RBFNODE *)malloc(n);
454 +                if (rbf == NULL)
455 +                        goto memerr;
456 +                memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
457 +                rbf->next = NULL; rbf->ejl = NULL;
458 +                return(rbf);
459 +        }
460 +        t /= full_dist;
461 + tryagain:
462 +        n = 0;                                  /* count migrating particles */
463 +        for (i = 0; i < mtx_nrows(mig); i++)
464 +            for (j = 0; j < mtx_ncols(mig); j++)
465 +                n += (mtx_coef(mig,i,j) > cthresh);
466 +                                                /* are we over our limit? */
467 +        if ((lobe_lim > 0) & (n > lobe_lim)) {
468 +                cthresh = cthresh*2. + 10.*FTINY;
469 +                goto tryagain;
470 +        }
471 + #ifdef DEBUG
472 +        fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
473 +                        mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
474 + #endif
475 +        rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
476 +        if (rbf == NULL)
477 +                goto memerr;
478 +        rbf->next = NULL; rbf->ejl = NULL;
479 +        VCOPY(rbf->invec, invec);
480 +        rbf->nrbf = n;
481 +        rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
482 +        n = 0;                                  /* advect RBF lobes */
483 +        for (i = 0; i < mtx_nrows(mig); i++) {
484 +            const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
485 +            const float         peak0 = rbf0i->peak;
486 +            const double        rad0 = R2ANG(rbf0i->crad);
487 +            FVECT               v0;
488 +            float               mv;
489 +            ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
490 +            for (j = 0; j < mtx_ncols(mig); j++)
491 +                if ((mv = mtx_coef(mig,i,j)) > cthresh) {
492 +                        const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
493 +                        double          rad2;
494 +                        FVECT           v;
495 +                        int             pos[2];
496 +                        rad2 = R2ANG(rbf1j->crad);
497 +                        rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
498 +                        rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
499 +                                                rad0*rad0/rad2;
500 +                        rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
501 +                        ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
502 +                        geodesic(v, v0, v, t, GEOD_REL);
503 +                        pos_from_vec(pos, v);
504 +                        rbf->rbfa[n].gx = pos[0];
505 +                        rbf->rbfa[n].gy = pos[1];
506 +                        ++n;
507 +                }
508 +        }
509 +        rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
510 +        return(rbf);
511 + memerr:
512 +        fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname);
513 +        exit(1);
514 +        return(NULL);   /* pro forma return */
515 + }
516 +
517   /* Clear our BSDF representation and free memory */
518   void
519   clear_bsdf_rep(void)
# Line 379 | Line 528 | clear_bsdf_rep(void)
528                  dsf_list = rbf->next;
529                  free(rbf);
530          }
531 +        bsdf_name[0] = '\0';
532 +        bsdf_manuf[0] = '\0';
533          inp_coverage = 0;
534          single_plane_incident = -1;
535          input_orient = output_orient = 0;
536          grid_res = GRIDRES;
537 +        bsdf_min = 0;
538 +        bsdf_spec_peak = 0;
539 +        bsdf_spec_rad = 0;
540   }
541  
542   /* Write our BSDF mesh interpolant out to the given binary stream */
# Line 393 | Line 547 | save_bsdf_rep(FILE *ofp)
547          MIGRATION       *mig;
548          int             i, n;
549                                          /* finish header */
550 +        if (bsdf_name[0])
551 +                fprintf(ofp, "NAME=%s\n", bsdf_name);
552 +        if (bsdf_manuf[0])
553 +                fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf);
554          fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
555          fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
556          fprintf(ofp, "GRIDRES=%d\n", grid_res);
557 +        fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
558 +        if ((bsdf_spec_peak > bsdf_min) & (bsdf_spec_rad > 0))
559 +                fprintf(ofp, "BSDFSPEC= %f %f\n", bsdf_spec_peak, bsdf_spec_rad);
560          fputformat(BSDFREP_FMT, ofp);
561          fputc('\n', ofp);
562                                          /* write each DSF */
# Line 448 | Line 609 | headline(char *s, void *p)
609   {
610          char    fmt[32];
611  
612 +        if (!strncmp(s, "NAME=", 5)) {
613 +                strcpy(bsdf_name, s+5);
614 +                bsdf_name[strlen(bsdf_name)-1] = '\0';
615 +        }
616 +        if (!strncmp(s, "MANUFACT=", 9)) {
617 +                strcpy(bsdf_manuf, s+9);
618 +                bsdf_manuf[strlen(bsdf_manuf)-1] = '\0';
619 +        }
620          if (!strncmp(s, "SYMMETRY=", 9)) {
621                  inp_coverage = atoi(s+9);
622                  single_plane_incident = !inp_coverage;
# Line 461 | Line 630 | headline(char *s, void *p)
630                  sscanf(s+8, "%d", &grid_res);
631                  return(0);
632          }
633 +        if (!strncmp(s, "BSDFMIN=", 8)) {
634 +                sscanf(s+8, "%lf", &bsdf_min);
635 +                return(0);
636 +        }
637 +        if (!strncmp(s, "BSDFSPEC=", 9)) {
638 +                sscanf(s+9, "%lf %lf", &bsdf_spec_peak, &bsdf_spec_rad);
639 +                return(0);
640 +        }
641          if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
642                  return(-1);
643          return(0);
# Line 477 | Line 654 | load_bsdf_rep(FILE *ifp)
654          clear_bsdf_rep();
655          if (ifp == NULL)
656                  return(0);
657 <        if (getheader(ifp, headline, NULL) < 0 || single_plane_incident < 0 |
658 <                        !input_orient | !output_orient) {
657 >        if (getheader(ifp, headline, NULL) < 0 || (single_plane_incident < 0) |
658 >                        !input_orient | !output_orient |
659 >                        (grid_res < 16) | (grid_res > 256)) {
660                  fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
661                                  progname);
662                  return(0);
663          }
664 <        rbfh.next = NULL;               /* read each DSF */
487 <        rbfh.ejl = NULL;
664 >        memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */
665          while ((rbfh.ord = getint(4, ifp)) >= 0) {
666                  RBFNODE         *newrbf;
667  
668                  rbfh.invec[0] = getflt(ifp);
669                  rbfh.invec[1] = getflt(ifp);
670                  rbfh.invec[2] = getflt(ifp);
671 +                if (normalize(rbfh.invec) == 0) {
672 +                        fprintf(stderr, "%s: zero incident vector\n", progname);
673 +                        return(0);
674 +                }
675                  rbfh.vtotal = getflt(ifp);
676                  rbfh.nrbf = getint(4, ifp);
677                  newrbf = (RBFNODE *)malloc(sizeof(RBFNODE) +
678                                          sizeof(RBFVAL)*(rbfh.nrbf-1));
679                  if (newrbf == NULL)
680                          goto memerr;
681 <                memcpy(newrbf, &rbfh, sizeof(RBFNODE));
681 >                *newrbf = rbfh;
682                  for (i = 0; i < rbfh.nrbf; i++) {
683                          newrbf->rbfa[i].peak = getflt(ifp);
684                          newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;

Diff Legend

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