--- ray/src/common/bsdf_t.c 2011/08/21 22:38:12 3.20 +++ ray/src/common/bsdf_t.c 2014/03/21 17:49:53 3.33 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf_t.c,v 3.20 2011/08/21 22:38:12 greg Exp $"; +static const char RCSid[] = "$Id: bsdf_t.c,v 3.33 2014/03/21 17:49:53 greg Exp $"; #endif /* * bsdf_t.c @@ -26,15 +26,18 @@ typedef int SDtreCallback(float val, const double *cmi /* reference width maximum (1.0) */ static const unsigned iwbits = sizeof(unsigned)*4; -static const unsigned iwmax = (1<<(sizeof(unsigned)*4))-1; +static const unsigned iwmax = 1<<(sizeof(unsigned)*4); /* maximum cumulative value */ static const unsigned cumlmax = ~0; /* constant z-vector */ static const FVECT zvec = {.0, .0, 1.}; + /* quantization value */ +static double quantum = 1./256.; /* Struct used for our distribution-building callback */ typedef struct { - int nic; /* number of input coordinates */ + short nic; /* number of input coordinates */ + short rev; /* reversing query */ unsigned alen; /* current array length */ unsigned nall; /* number of allocated entries */ unsigned wmin; /* minimum square size so far */ @@ -281,7 +284,7 @@ SDdotravTre(const SDNode *st, const double *pos, int c unsigned skipmask = 0; csiz *= .5; for (i = st->ndim; i--; ) - if (1<ndim; n--; ) { if (n & 1<ndim; n--; ) { if (1<sidef) { /* whose side are you on? */ - case SD_UFRONT: + case SD_FREFL: if ((outVec[2] < 0) | (inVec[2] < 0)) return -1.; break; - case SD_UBACK: + case SD_BREFL: if ((outVec[2] > 0) | (inVec[2] > 0)) return -1.; break; - case SD_XMIT: - if ((outVec[2] > 0) == (inVec[2] > 0)) + case SD_FXMIT: + if (outVec[2] > 0) { + if (inVec[2] > 0) + return -1.; + vtmp = outVec; outVec = inVec; inVec = vtmp; + } else if (inVec[2] < 0) return -1.; break; + case SD_BXMIT: + if (inVec[2] > 0) { + if (outVec[2] > 0) + return -1.; + vtmp = outVec; outVec = inVec; inVec = vtmp; + } else if (outVec[2] < 0) + return -1.; + break; default: return -1.; } /* convert vector coordinates */ if (sdt->st->ndim == 3) { spinvector(rOutVec, outVec, zvec, -atan2(-inVec[1],-inVec[0])); - gridPos[0] = .5 - .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]); + gridPos[0] = (.5-FTINY) - + .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]); SDdisk2square(gridPos+1, rOutVec[0], rOutVec[1]); } else if (sdt->st->ndim == 4) { SDdisk2square(gridPos, -inVec[0], -inVec[1]); @@ -483,9 +501,16 @@ build_scaffold(float val, const double *cmin, double c { SDdistScaffold *sp = (SDdistScaffold *)cptr; int wid = csiz*(double)iwmax + .5; + double revcmin[2]; bitmask_t bmin[2], bmax[2]; - cmin += sp->nic; /* skip to output coords */ + if (sp->rev) { /* need to reverse sense? */ + revcmin[0] = 1. - cmin[0] - csiz; + revcmin[1] = 1. - cmin[1] - csiz; + cmin = revcmin; + } else { + cmin += sp->nic; /* else skip to output coords */ + } if (wid < sp->wmin) /* new minimum width? */ sp->wmin = wid; if (wid > sp->wmax) /* new maximum? */ @@ -531,9 +556,11 @@ sscmp(const void *p1, const void *p2) /* Create a new cumulative distribution for the given input direction */ static SDTreCDst * -make_cdist(const SDTre *sdt, const double *pos) +make_cdist(const SDTre *sdt, const double *invec, int rev) { SDdistScaffold myScaffold; + double pos[4]; + int cmask; SDTreCDst *cd; struct outdir_s *sp; double scale, cursum; @@ -542,14 +569,20 @@ make_cdist(const SDTre *sdt, const double *pos) myScaffold.wmin = iwmax; myScaffold.wmax = 0; myScaffold.nic = sdt->st->ndim - 2; + myScaffold.rev = rev; myScaffold.alen = 0; myScaffold.nall = 512; myScaffold.darr = (struct outdir_s *)malloc(sizeof(struct outdir_s) * myScaffold.nall); if (myScaffold.darr == NULL) return NULL; + /* set up traversal */ + cmask = (1<st, pos, (1<st, pos, cmask, &build_scaffold, &myScaffold) < 0) { free(myScaffold.darr); return NULL; @@ -572,7 +605,7 @@ make_cdist(const SDTre *sdt, const double *pos) /* record input range */ scale = myScaffold.wmin / (double)iwmax; for (i = myScaffold.nic; i--; ) { - cd->clim[i][0] = floor(pos[i]/scale) * scale; + cd->clim[i][0] = floor(pos[i+2*rev]/scale) * scale; cd->clim[i][1] = cd->clim[i][0] + scale; } if (cd->isodist) { /* avoid issue in SDqueryTreProjSA() */ @@ -581,7 +614,10 @@ make_cdist(const SDTre *sdt, const double *pos) } cd->max_psa = myScaffold.wmax / (double)iwmax; cd->max_psa *= cd->max_psa * M_PI; - cd->sidef = sdt->sidef; + if (rev) + cd->sidef = (sdt->sidef==SD_BXMIT) ? SD_FXMIT : SD_BXMIT; + else + cd->sidef = sdt->sidef; cd->cTotal = 1e-20; /* compute directional total */ sp = myScaffold.darr; for (i = myScaffold.alen; i--; sp++) @@ -609,20 +645,52 @@ SDgetTreCDist(const FVECT inVec, SDComponent *sdc) const SDTre *sdt; double inCoord[2]; int i; + int mode; SDTreCDst *cd, *cdlast; /* check arguments */ if ((inVec == NULL) | (sdc == NULL) || (sdt = (SDTre *)sdc->dist) == NULL) return NULL; - if (sdt->st->ndim == 3) /* isotropic BSDF? */ - inCoord[0] = .5 - .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]); - else if (sdt->st->ndim == 4) - SDdisk2square(inCoord, -inVec[0], -inVec[1]); - else + switch (mode = sdt->sidef) { /* check direction */ + case SD_FREFL: + if (inVec[2] < 0) + return NULL; + break; + case SD_BREFL: + if (inVec[2] > 0) + return NULL; + break; + case SD_FXMIT: + if (inVec[2] < 0) + mode = SD_BXMIT; + break; + case SD_BXMIT: + if (inVec[2] > 0) + mode = SD_FXMIT; + break; + default: + return NULL; + } + if (sdt->st->ndim == 3) { /* isotropic BSDF? */ + if (mode != sdt->sidef) /* XXX unhandled reciprocity */ + return &SDemptyCD; + inCoord[0] = (.5-FTINY) - + .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]); + } else if (sdt->st->ndim == 4) { + if (mode != sdt->sidef) /* use reciprocity? */ + SDdisk2square(inCoord, inVec[0], inVec[1]); + else + SDdisk2square(inCoord, -inVec[0], -inVec[1]); + } else return NULL; /* should be internal error */ + /* quantize to avoid f.p. errors */ + for (i = sdt->st->ndim - 2; i--; ) + inCoord[i] = floor(inCoord[i]/quantum)*quantum + .5*quantum; cdlast = NULL; /* check for direction in cache list */ for (cd = (SDTreCDst *)sdc->cdList; cd != NULL; cdlast = cd, cd = cd->next) { + if (cd->sidef != mode) + continue; for (i = sdt->st->ndim - 2; i--; ) if ((cd->clim[i][0] > inCoord[i]) | (inCoord[i] >= cd->clim[i][1])) @@ -631,7 +699,7 @@ SDgetTreCDist(const FVECT inVec, SDComponent *sdc) break; /* means we have a match */ } if (cd == NULL) /* need to create new entry? */ - cdlast = cd = make_cdist(sdt, inCoord); + cdlast = cd = make_cdist(sdt, inCoord, mode != sdt->sidef); if (cdlast != NULL) { /* move entry to head of cache list */ cdlast->next = cd->next; cd->next = (SDTreCDst *)sdc->cdList; @@ -663,10 +731,12 @@ SDqueryTreProjSA(double *psa, const FVECT v1, const RR } else { const SDTreCDst *cd = (const SDTreCDst *)SDgetTreCDist(v1, sdc); if (cd == NULL) - return SDEmemory; - myPSA[0] = M_PI * (cd->clim[0][1] - cd->clim[0][0]) * - (cd->clim[1][1] - cd->clim[1][0]); - myPSA[1] = cd->max_psa; + myPSA[0] = myPSA[1] = 0; + else { + myPSA[0] = M_PI * (cd->clim[0][1] - cd->clim[0][0]) * + (cd->clim[1][1] - cd->clim[1][0]); + myPSA[1] = cd->max_psa; + } } switch (qflags) { /* record based on flag settings */ case SDqueryVal: @@ -681,7 +751,7 @@ SDqueryTreProjSA(double *psa, const FVECT v1, const RR psa[1] = myPSA[1]; /* fall through */ case SDqueryMin: - if (myPSA[0] < psa[0]) + if ((myPSA[0] > 0) & (myPSA[0] < psa[0])) psa[0] = myPSA[0]; break; } @@ -702,10 +772,12 @@ SDsampTreCDist(FVECT ioVec, double randX, const SDCDst /* check arguments */ if ((ioVec == NULL) | (cd == NULL)) return SDEargument; + if (!cd->sidef) + return SDEnone; /* XXX should never happen */ if (ioVec[2] > 0) { - if (!(cd->sidef & SD_UFRONT)) + if ((cd->sidef != SD_FREFL) & (cd->sidef != SD_FXMIT)) return SDEargument; - } else if (!(cd->sidef & SD_UBACK)) + } else if ((cd->sidef != SD_BREFL) & (cd->sidef != SD_BXMIT)) return SDEargument; /* binary search to find position */ ilower = 0; iupper = cd->calen; @@ -728,10 +800,9 @@ SDsampTreCDist(FVECT ioVec, double randX, const SDCDst SDsquare2disk(gpos, gpos[0], gpos[1]); /* compute Z-coordinate */ gpos[2] = 1. - gpos[0]*gpos[0] - gpos[1]*gpos[1]; - if (gpos[2] > 0) /* paranoia, I hope */ - gpos[2] = sqrt(gpos[2]); + gpos[2] = sqrt(gpos[2]*(gpos[2]>0)); /* emit from back? */ - if (ioVec[2] > 0 ^ cd->sidef != SD_XMIT) + if ((cd->sidef == SD_BREFL) | (cd->sidef == SD_FXMIT)) gpos[2] = -gpos[2]; if (cd->isodist) { /* rotate isotropic result */ rotangle = atan2(-ioVec[1],-ioVec[0]); @@ -778,7 +849,8 @@ load_values(char **spp, float *va, int n) char *svnext; while (n-- > 0 && (svnext = fskip(*spp)) != NULL) { - *v++ = atof(*spp); + if ((*v++ = atof(*spp)) < 0) + v[-1] = 0; *spp = svnext; eat_token(spp, ','); } @@ -841,6 +913,8 @@ get_extrema(SDSpectralDF *df) double stepWidth, dhemi, bmin[4], bmax[4]; stepWidth = SDsmallestLeaf(st); + if (quantum > stepWidth) /* adjust quantization factor */ + quantum = stepWidth; df->minProjSA = M_PI*stepWidth*stepWidth; if (stepWidth < .03125) stepWidth = .03125; /* 1/32 resolution good enough */ @@ -887,20 +961,26 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim) /* * Remember that front and back are reversed from WINDOW 6 orientations */ - if (!strcasecmp(sdata, "Transmission")) { + if (!strcasecmp(sdata, "Transmission Front")) { + if (sd->tb != NULL) + SDfreeSpectralDF(sd->tb); + if ((sd->tb = SDnewSpectralDF(1)) == NULL) + return SDEmemory; + df = sd->tb; + } else if (!strcasecmp(sdata, "Transmission Back")) { if (sd->tf != NULL) SDfreeSpectralDF(sd->tf); if ((sd->tf = SDnewSpectralDF(1)) == NULL) return SDEmemory; df = sd->tf; } else if (!strcasecmp(sdata, "Reflection Front")) { - if (sd->rb != NULL) /* note back-front reversal */ + if (sd->rb != NULL) SDfreeSpectralDF(sd->rb); if ((sd->rb = SDnewSpectralDF(1)) == NULL) return SDEmemory; df = sd->rb; } else if (!strcasecmp(sdata, "Reflection Back")) { - if (sd->rf != NULL) /* note front-back reversal */ + if (sd->rf != NULL) SDfreeSpectralDF(sd->rf); if ((sd->rf = SDnewSpectralDF(1)) == NULL) return SDEmemory; @@ -920,11 +1000,13 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim) if (sdt == NULL) return SDEmemory; if (df == sd->rf) - sdt->sidef = SD_UFRONT; + sdt->sidef = SD_FREFL; else if (df == sd->rb) - sdt->sidef = SD_UBACK; - else - sdt->sidef = SD_XMIT; + sdt->sidef = SD_BREFL; + else if (df == sd->tf) + sdt->sidef = SD_FXMIT; + else /* df == sd->tb */ + sdt->sidef = SD_BXMIT; sdt->st = NULL; df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */ df->comp[0].dist = sdt; @@ -1086,7 +1168,10 @@ SDloadTre(SDData *sd, ezxml_t wtl) /* separate diffuse components */ extract_diffuse(&sd->rLambFront, sd->rf); extract_diffuse(&sd->rLambBack, sd->rb); - extract_diffuse(&sd->tLamb, sd->tf); + if (sd->tf != NULL) + extract_diffuse(&sd->tLamb, sd->tf); + if (sd->tb != NULL) + extract_diffuse(&sd->tLamb, sd->tb); /* return success */ return SDEnone; }