--- ray/src/common/bsdf_t.c 2011/05/01 16:34:37 3.12 +++ ray/src/common/bsdf_t.c 2018/05/10 22:55:35 3.47 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf_t.c,v 3.12 2011/05/01 16:34:37 greg Exp $"; +static const char RCSid[] = "$Id: bsdf_t.c,v 3.47 2018/05/10 22:55:35 greg Exp $"; #endif /* * bsdf_t.c @@ -10,6 +10,7 @@ static const char RCSid[] = "$Id: bsdf_t.c,v 3.12 2011 * */ +#define _USE_MATH_DEFINES #include "rtio.h" #include #include @@ -20,18 +21,29 @@ static const char RCSid[] = "$Id: bsdf_t.c,v 3.12 2011 #include "hilbert.h" /* Callback function type for SDtraverseTre() */ -typedef int SDtreCallback(float val, const double *cmin, - double csiz, void *cptr); - +typedef int SDtreCallback(float val, const double *cmin, double csiz, + void *cptr); /* 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.; + /* our RGB primaries */ +static C_COLOR tt_RGB_prim[3]; +static float tt_RGB_coef[3]; +static const double czero[SD_MAXDIM]; + +enum {tt_Y, tt_u, tt_v}; /* tree components (tt_Y==0) */ + /* 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 */ @@ -91,7 +103,7 @@ SDfreeTre(SDNode *st) return; for (n = (st->log2GR < 0) << st->ndim; n--; ) SDfreeTre(st->u.t[n]); - free((void *)st); + free(st); } /* Free a variable-resolution BSDF */ @@ -102,7 +114,9 @@ SDFreeBTre(void *p) if (sdt == NULL) return; - SDfreeTre(sdt->st); + SDfreeTre(sdt->stc[tt_Y]); + SDfreeTre(sdt->stc[tt_u]); + SDfreeTre(sdt->stc[tt_v]); free(sdt); } @@ -125,13 +139,13 @@ fill_grid_branch(float *dptr, const float *sptr, int n static float * grid_branch_start(SDNode *st, int n) { - unsigned skipsiz = 1 << st->log2GR; + unsigned skipsiz = 1 << (st->log2GR - 1); float *vptr = st->u.v; int i; - for (i = 0; i < st->ndim; skipsiz <<= st->log2GR) - if (1<> 1; + for (i = st->ndim; i--; skipsiz <<= st->log2GR) + if (1< 0; ) + if ((tmin[i] < .0) | (tmin[i] >= 1.-FTINY)) + break; + if ((i >= 0) | (nd <= 0) | (tsiz <= FTINY) | (tsiz > 1.+FTINY) | + (sroot != NULL && sroot->ndim != nd)) { + SDfreeTre(sroot); + return NULL; + } + if (tsiz >= 1.-FTINY) { /* special case when tree is a leaf */ + SDfreeTre(sroot); + if ((sroot = SDnewNode(nd, 0)) != NULL) + sroot->u.v[0] = val; + return sroot; + } + /* make sure we have branching root */ + if (sroot != NULL && sroot->log2GR >= 0) { + SDfreeTre(sroot); sroot = NULL; + } + if (sroot == NULL && (sroot = SDnewNode(nd, -1)) == NULL) + return NULL; + st = sroot; /* climb/grow tree */ + memset(ctrk, 0, sizeof(ctrk)); + for ( ; ; ) { + csiz *= .5; /* find appropriate branch */ + n = 0; + for (i = nd; i--; ) + if (ctrk[i]+csiz <= tmin[i]+FTINY) { + ctrk[i] += csiz; + n |= 1 << i; + } + /* reached desired voxel? */ + if (csiz <= tsiz+FTINY) { + SDfreeTre(st->u.t[n]); + st = st->u.t[n] = SDnewNode(nd, 0); + break; + } + /* else grow tree as needed */ + if (st->u.t[n] != NULL && st->u.t[n]->log2GR >= 0) { + SDfreeTre(st->u.t[n]); st->u.t[n] = NULL; + } + if (st->u.t[n] == NULL) + st->u.t[n] = SDnewNode(nd, -1); + if ((st = st->u.t[n]) == NULL) + break; + } + if (st == NULL) { + SDfreeTre(sroot); + return NULL; + } + st->u.v[0] = val; /* assign leaf and return root */ + return sroot; +} + /* Find smallest leaf in tree */ static double SDsmallestLeaf(const SDNode *st) @@ -190,13 +266,15 @@ SDiterSum(const float *va, int nd, int shft, const int const unsigned skipsiz = 1 << --nd*shft; double sum = .0; int i; - + + va += *imin * skipsiz; + if (skipsiz == 1) for (i = *imin; i < *imax; i++) - sum += va[i]; + sum += *va++; else - for (i = *imin; i < *imax; i++) - sum += SDiterSum(va + i*skipsiz, nd, shft, imin+1, imax+1); + for (i = *imin; i < *imax; i++, va += skipsiz) + sum += SDiterSum(va, nd, shft, imin+1, imax+1); return sum; } @@ -204,7 +282,6 @@ SDiterSum(const float *va, int nd, int shft, const int static double SDavgTreBox(const SDNode *st, const double *bmin, const double *bmax) { - int imin[SD_MAXDIM], imax[SD_MAXDIM]; unsigned n; int i; @@ -214,7 +291,7 @@ SDavgTreBox(const SDNode *st, const double *bmin, cons for (i = st->ndim; i--; ) { if (bmin[i] >= 1.) return .0; - if (bmax[i] <= .0) + if (bmax[i] <= 0) return .0; if (bmin[i] >= bmax[i]) return .0; @@ -222,7 +299,6 @@ SDavgTreBox(const SDNode *st, const double *bmin, cons if (st->log2GR < 0) { /* iterate on subtree */ double sum = .0, wsum = 1e-20; double sbmin[SD_MAXDIM], sbmax[SD_MAXDIM], w; - for (n = 1 << st->ndim; n--; ) { w = 1.; for (i = st->ndim; i--; ) { @@ -234,6 +310,10 @@ SDavgTreBox(const SDNode *st, const double *bmin, cons } if (sbmin[i] < .0) sbmin[i] = .0; if (sbmax[i] > 1.) sbmax[i] = 1.; + if (sbmin[i] >= sbmax[i]) { + w = .0; + break; + } w *= sbmax[i] - sbmin[i]; } if (w > 1e-10) { @@ -242,19 +322,22 @@ SDavgTreBox(const SDNode *st, const double *bmin, cons } } return sum / wsum; + } else { /* iterate over leaves */ + int imin[SD_MAXDIM], imax[SD_MAXDIM]; + + n = 1; + for (i = st->ndim; i--; ) { + imin[i] = (bmin[i] <= 0) ? 0 : + (int)((1 << st->log2GR)*bmin[i]); + imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR) : + (int)((1 << st->log2GR)*bmax[i] + .999999); + n *= imax[i] - imin[i]; + } + if (n) + return SDiterSum(st->u.v, st->ndim, + st->log2GR, imin, imax) / (double)n; } - n = 1; /* iterate over leaves */ - for (i = st->ndim; i--; ) { - imin[i] = (bmin[i] <= 0) ? 0 - : (int)((1 << st->log2GR)*bmin[i]); - imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR) - : (int)((1 << st->log2GR)*bmax[i] + .999999); - n *= imax[i] - imin[i]; - } - if (!n) - return .0; - - return SDiterSum(st->u.v, st->ndim, st->log2GR, imin, imax) / (double)n; + return .0; } /* Recursive call for SDtraverseTre() */ @@ -266,21 +349,26 @@ SDdotravTre(const SDNode *st, const double *pos, int c int rv, rval = 0; double bmin[SD_MAXDIM]; int i, n; + /* paranoia */ + if (st == NULL) + return 0; /* in branches? */ if (st->log2GR < 0) { unsigned skipmask = 0; - csiz *= .5; for (i = st->ndim; i--; ) - if (1<ndim; n--; ) + for (n = 1 << st->ndim; n--; ) { if (n & 1<ndim; n--; ) + for (n = 1 << st->ndim; n--; ) { if (!(n & 1<ndim; n--; ) { if (1<log2GR; } - /* fill in unused dimensions */ - for (i = SD_MAXDIM; i-- > st->ndim; ) { - clim[i][0] = 0; clim[i][1] = 1; - } #if (SD_MAXDIM == 4) bmin[0] = cmin[0] + csiz*clim[0][0]; for (cpos[0] = clim[0][0]; cpos[0] < clim[0][1]; cpos[0]++) { bmin[1] = cmin[1] + csiz*clim[1][0]; for (cpos[1] = clim[1][0]; cpos[1] < clim[1][1]; cpos[1]++) { bmin[2] = cmin[2] + csiz*clim[2][0]; - for (cpos[2] = clim[2][0]; cpos[2] < clim[2][1]; cpos[2]++) { - bmin[3] = cmin[3] + csiz*(cpos[3] = clim[3][0]); + if (st->ndim == 3) { + cpos[2] = clim[2][0]; n = cpos[0]; - for (i = 1; i < st->ndim; i++) + for (i = 1; i < 3; i++) n = (n << st->log2GR) + cpos[i]; - for ( ; cpos[3] < clim[3][1]; cpos[3]++) { + for ( ; cpos[2] < clim[2][1]; cpos[2]++) { rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr); if (rv < 0) return rv; - bmin[3] += csiz; + bmin[2] += csiz; } - bmin[2] += csiz; + } else { + for (cpos[2] = clim[2][0]; cpos[2] < clim[2][1]; cpos[2]++) { + bmin[3] = cmin[3] + csiz*(cpos[3] = clim[3][0]); + n = cpos[0]; + for (i = 1; i < 4; i++) + n = (n << st->log2GR) + cpos[i]; + for ( ; cpos[3] < clim[3][1]; cpos[3]++) { + rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr); + if (rv < 0) + return rv; + bmin[3] += csiz; + } + bmin[2] += csiz; + } } bmin[1] += csiz; } @@ -353,7 +450,6 @@ static int SDtraverseTre(const SDNode *st, const double *pos, int cmask, SDtreCallback *cf, void *cptr) { - static double czero[SD_MAXDIM]; int i; /* check arguments */ if ((st == NULL) | (cf == NULL)) @@ -378,7 +474,7 @@ SDlookupTre(const SDNode *st, const double *pos, doubl hcube[i] = .0; } /* climb the tree */ - while (st->log2GR < 0) { + while (st != NULL && st->log2GR < 0) { n = 0; /* move to appropriate branch */ if (hcube) hcube[st->ndim] *= .5; for (i = st->ndim; i--; ) { @@ -391,6 +487,8 @@ SDlookupTre(const SDNode *st, const double *pos, doubl st = st->u.t[n]; /* avoids tail recursion */ pos = spos; } + if (st == NULL) /* should never happen? */ + return .0; if (st->log2GR == 0) /* short cut */ return st->u.v[0]; n = t = 0; /* find grid array index */ @@ -406,42 +504,84 @@ SDlookupTre(const SDNode *st, const double *pos, doubl return st->u.v[n]; /* no interpolation */ } +/* Convert CIE (Y,u',v') color to our RGB */ +static void +SDyuv2rgb(double yval, double uprime, double vprime, float rgb[3]) +{ + const double dfact = 1./(6.*uprime - 16.*vprime + 12.); + C_COLOR cxy; + + c_cset(&cxy, 9.*uprime*dfact, 4.*vprime*dfact); + c_toSharpRGB(&cxy, yval, rgb); +} + /* Query BSDF value and sample hypercube for the given vectors */ -static float -SDqueryTre(const SDTre *sdt, const FVECT outVec, const FVECT inVec, double *hc) +static int +SDqueryTre(const SDTre *sdt, float *coef, + const FVECT outVec, const FVECT inVec, double *hc) { - static const FVECT zvec = {.0, .0, 1.}; - FVECT rOutVec; - double gridPos[4]; + const RREAL *vtmp; + float yval; + FVECT rOutVec; + double gridPos[4]; + if (sdt->stc[tt_Y] == NULL) /* paranoia, I hope */ + return 0; + switch (sdt->sidef) { /* whose side are you on? */ - case SD_UFRONT: + case SD_FREFL: if ((outVec[2] < 0) | (inVec[2] < 0)) - return -1.; + return 0; break; - case SD_UBACK: + case SD_BREFL: if ((outVec[2] > 0) | (inVec[2] > 0)) - return -1.; + return 0; break; - case SD_XMIT: - if ((outVec[2] > 0) == (inVec[2] > 0)) - return -1.; + case SD_FXMIT: + if (outVec[2] > 0) { + if (inVec[2] > 0) + return 0; + vtmp = outVec; outVec = inVec; inVec = vtmp; + } else if (inVec[2] < 0) + return 0; break; + case SD_BXMIT: + if (inVec[2] > 0) { + if (outVec[2] > 0) + return 0; + vtmp = outVec; outVec = inVec; inVec = vtmp; + } else if (outVec[2] < 0) + return 0; + break; default: - return -1.; + return 0; } /* 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]); + if (sdt->stc[tt_Y]->ndim == 3) { + spinvector(rOutVec, outVec, zvec, -atan2(-inVec[1],-inVec[0])); + 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) { + } else if (sdt->stc[tt_Y]->ndim == 4) { SDdisk2square(gridPos, -inVec[0], -inVec[1]); SDdisk2square(gridPos+2, outVec[0], outVec[1]); } else - return -1.; /* should be internal error */ - - return SDlookupTre(sdt->st, gridPos, hc); + return 0; /* should be internal error */ + /* get BSDF value */ + yval = SDlookupTre(sdt->stc[tt_Y], gridPos, hc); + if (coef == NULL) /* just getting hypercube? */ + return 1; + if (sdt->stc[tt_u] == NULL || sdt->stc[tt_v] == NULL) { + *coef = yval; + return 1; /* no color */ + } + /* else decode color */ + SDyuv2rgb(yval, SDlookupTre(sdt->stc[tt_u], gridPos, NULL), + SDlookupTre(sdt->stc[tt_v], gridPos, NULL), coef); + coef[0] *= tt_RGB_coef[0]; + coef[1] *= tt_RGB_coef[1]; + coef[2] *= tt_RGB_coef[2]; + return 3; } /* Compute non-diffuse component for variable-resolution BSDF */ @@ -454,8 +594,7 @@ SDgetTreBSDF(float coef[SDmaxCh], const FVECT outVec, || sdc->dist == NULL) return 0; /* get nearest BSDF value */ - coef[0] = SDqueryTre((SDTre *)sdc->dist, outVec, inVec, NULL); - return (coef[0] >= 0); /* monochromatic for now */ + return SDqueryTre((SDTre *)sdc->dist, coef, outVec, inVec, NULL); } /* Callback to build cumulative distribution using SDtraverseTre() */ @@ -464,16 +603,23 @@ 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? */ sp->wmax = wid; if (sp->alen >= sp->nall) { /* need more space? */ struct outdir_s *ndarr; - sp->nall += 1024; + sp->nall = (int)(1.5*sp->nall) + 256; ndarr = (struct outdir_s *)realloc(sp->darr, sizeof(struct outdir_s)*sp->nall); if (ndarr == NULL) { @@ -512,9 +658,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; @@ -522,16 +670,22 @@ make_cdist(const SDTre *sdt, const double *pos) /* initialize scaffold */ myScaffold.wmin = iwmax; myScaffold.wmax = 0; - myScaffold.nic = sdt->st->ndim - 2; + myScaffold.nic = sdt->stc[tt_Y]->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<stc[tt_Y], pos, cmask, + build_scaffold, &myScaffold) < 0) { free(myScaffold.darr); return NULL; } @@ -545,19 +699,27 @@ make_cdist(const SDTre *sdt, const double *pos) free(myScaffold.darr); return NULL; } + cd->isodist = (myScaffold.nic == 1); /* sort the distribution */ qsort(myScaffold.darr, cd->calen = myScaffold.alen, - sizeof(struct outdir_s), &sscmp); + sizeof(struct outdir_s), sscmp); /* 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() */ + cd->clim[1][0] = cd->clim[0][0]; + cd->clim[1][1] = cd->clim[0][1]; + } 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++) @@ -584,23 +746,55 @@ SDgetTreCDist(const FVECT inVec, SDComponent *sdc) { const SDTre *sdt; double inCoord[2]; - int vflags; 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->stc[tt_Y]->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->stc[tt_Y]->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->stc[tt_Y]->ndim - 2; i--; ) + inCoord[i] = floor(inCoord[i]/quantum)*quantum + .5*quantum; cdlast = NULL; /* check for direction in cache list */ + /* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */ for (cd = (SDTreCDst *)sdc->cdList; cd != NULL; - cdlast = cd, cd = (SDTreCDst *)cd->next) { - for (i = sdt->st->ndim - 2; i--; ) + cdlast = cd, cd = cd->next) { + if (cd->sidef != mode) + continue; + for (i = sdt->stc[tt_Y]->ndim - 2; i--; ) if ((cd->clim[i][0] > inCoord[i]) | (inCoord[i] >= cd->clim[i][1])) break; @@ -608,12 +802,13 @@ 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 = sdc->cdList; + cd->next = (SDTreCDst *)sdc->cdList; sdc->cdList = (SDCDst *)cd; } + /* END MUTEX LOCK */ return (SDCDst *)cd; /* ready to go */ } @@ -630,20 +825,22 @@ SDqueryTreProjSA(double *psa, const FVECT v1, const RR /* get projected solid angle(s) */ if (v2 != NULL) { const SDTre *sdt = (SDTre *)sdc->dist; - double hcube[SD_MAXDIM]; - if (SDqueryTre(sdt, v1, v2, hcube) < 0) { + double hcube[SD_MAXDIM+1]; + if (!SDqueryTre(sdt, NULL, v1, v2, hcube)) { strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA"); return SDEinternal; } - myPSA[0] = hcube[sdt->st->ndim]; + myPSA[0] = hcube[sdt->stc[tt_Y]->ndim]; myPSA[1] = myPSA[0] *= myPSA[0] * M_PI; } 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: @@ -658,7 +855,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; } @@ -674,20 +871,22 @@ SDsampTreCDist(FVECT ioVec, double randX, const SDCDst const SDTreCDst *cd = (const SDTreCDst *)cdp; const unsigned target = randX*cumlmax; bitmask_t hndx, hcoord[2]; - double gpos[3]; + double gpos[3], rotangle; int i, iupper, ilower; /* 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; while ((i = (iupper + ilower) >> 1) != ilower) - if ((long)target >= (long)cd->carr[i].cuml) + if (target >= cd->carr[i].cuml) ilower = i; else iupper = i; @@ -705,12 +904,15 @@ 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]; - VCOPY(ioVec, gpos); + if (cd->isodist) { /* rotate isotropic sample */ + rotangle = atan2(-ioVec[1],-ioVec[0]); + spinvector(ioVec, gpos, zvec, rotangle); + } else + VCOPY(ioVec, gpos); return SDEnone; } @@ -724,7 +926,7 @@ next_token(char **spp) } /* Advance pointer past matching token (or any token if c==0) */ -#define eat_token(spp,c) (next_token(spp)==(c) ^ !(c) ? *(*(spp))++ : 0) +#define eat_token(spp,c) ((next_token(spp)==(c)) ^ !(c) ? *(*(spp))++ : 0) /* Count words from this point in string to '}' */ static int @@ -750,7 +952,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, ','); } @@ -780,7 +983,7 @@ load_tree_data(char **spp, int nd) } else { /* else load value grid */ int bsiz; n = count_values(*spp); /* see how big the grid is */ - for (bsiz = 0; bsiz < 8*sizeof(size_t)-1; bsiz += nd) + for (bsiz = 0; bsiz < 8*sizeof(size_t); bsiz += nd) if (1<= 8*sizeof(size_t)) { @@ -809,10 +1012,12 @@ load_tree_data(char **spp, int nd) static SDError get_extrema(SDSpectralDF *df) { - SDNode *st = (*(SDTre *)df->comp[0].dist).st; + SDNode *st = (*(SDTre *)df->comp[0].dist).stc[tt_Y]; 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 */ @@ -847,12 +1052,11 @@ get_extrema(SDSpectralDF *df) /* Load BSDF distribution for this wavelength */ static SDError -load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim) +load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int ndim) { SDSpectralDF *df; SDTre *sdt; char *sdata; - int i; /* allocate BSDF component */ sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection")); if (!sdata) @@ -860,27 +1064,24 @@ 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 (sd->tf != NULL) - SDfreeSpectralDF(sd->tf); - if ((sd->tf = SDnewSpectralDF(1)) == NULL) + if (!strcasecmp(sdata, "Transmission Front")) { + if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(1)) == NULL) return SDEmemory; + df = sd->tb; + } else if (!strcasecmp(sdata, "Transmission Back")) { + if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(1)) == NULL) + return SDEmemory; df = sd->tf; } else if (!strcasecmp(sdata, "Reflection Front")) { - if (sd->rb != NULL) /* note back-front reversal */ - SDfreeSpectralDF(sd->rb); - if ((sd->rb = SDnewSpectralDF(1)) == NULL) + if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(1)) == NULL) return SDEmemory; df = sd->rb; } else if (!strcasecmp(sdata, "Reflection Back")) { - if (sd->rf != NULL) /* note front-back reversal */ - SDfreeSpectralDF(sd->rf); - if ((sd->rf = SDnewSpectralDF(1)) == NULL) + if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(1)) == NULL) return SDEmemory; df = sd->rf; } else return SDEnone; - /* XXX should also check "ScatteringDataType" for consistency? */ /* get angle bases */ sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis")); if (!sdata || strcasecmp(sdata, "LBNL/Shirley-Chiu")) { @@ -888,20 +1089,28 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim) !sdata ? "Missing" : "Unsupported", sd->name); return !sdata ? SDEformat : SDEsupport; } - /* allocate BSDF tree */ - sdt = (SDTre *)malloc(sizeof(SDTre)); - if (sdt == NULL) - return SDEmemory; - if (df == sd->rf) - sdt->sidef = SD_UFRONT; - else if (df == sd->rb) - sdt->sidef = SD_UBACK; - else - sdt->sidef = SD_XMIT; - sdt->st = NULL; - df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */ - df->comp[0].dist = sdt; - df->comp[0].func = &SDhandleTre; + if (df->comp[0].dist == NULL) { /* need to allocate BSDF tree? */ + sdt = (SDTre *)malloc(sizeof(SDTre)); + if (sdt == NULL) + return SDEmemory; + if (df == sd->rf) + sdt->sidef = SD_FREFL; + else if (df == sd->rb) + sdt->sidef = SD_BREFL; + else if (df == sd->tf) + sdt->sidef = SD_FXMIT; + else /* df == sd->tb */ + sdt->sidef = SD_BXMIT; + sdt->stc[tt_Y] = sdt->stc[tt_u] = sdt->stc[tt_v] = NULL; + df->comp[0].dist = sdt; + df->comp[0].func = &SDhandleTre; + } else { + sdt = (SDTre *)df->comp[0].dist; + if (sdt->stc[ct] != NULL) { + SDfreeTre(sdt->stc[ct]); + sdt->stc[ct] = NULL; + } + } /* read BSDF data */ sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData")); if (!sdata || !next_token(&sdata)) { @@ -909,8 +1118,8 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim) sd->name); return SDEformat; } - sdt->st = load_tree_data(&sdata, ndim); - if (sdt->st == NULL) + sdt->stc[ct] = load_tree_data(&sdata, ndim); + if (sdt->stc[ct] == NULL) return SDEformat; if (next_token(&sdata)) { /* check for unconsumed characters */ sprintf(SDerrorDetail, @@ -919,10 +1128,11 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int ndim) return SDEformat; } /* flatten branches where possible */ - sdt->st = SDsimplifyTre(sdt->st); - if (sdt->st == NULL) + sdt->stc[ct] = SDsimplifyTre(sdt->stc[ct]); + if (sdt->stc[ct] == NULL) return SDEinternal; - return get_extrema(df); /* compute global quantities */ + /* compute global quantities for Y */ + return (ct == tt_Y) ? get_extrema(df) : SDEnone; } /* Find minimum value in tree */ @@ -957,21 +1167,23 @@ SDsubtractTreVal(SDNode *st, float val) SDsubtractTreVal(st->u.t[n], val); } else { for (n = 1<<(st->ndim*st->log2GR); n--; ) - st->u.v[n] -= val; + if ((st->u.v[n] -= val) < 0) + st->u.v[n] = .0f; } } -/* Subtract minimum value from BSDF */ +/* Subtract minimum Y value from BSDF */ static double -subtract_min(SDNode *st) +subtract_min_Y(SDNode *st) { - float vmin; + const float vmaxmin = 1.5/M_PI; + float vmin; /* be sure to skip unused portion */ if (st->ndim == 3) { int n; - vmin = 1./M_PI; + vmin = vmaxmin; if (st->log2GR < 0) { - for (n = 0; n < 4; n++) { + for (n = 0; n < 8; n += 2) { float v = SDgetTreMin(st->u.t[n]); if (v < vmin) vmin = v; @@ -985,36 +1197,172 @@ subtract_min(SDNode *st) } else /* anisotropic covers entire tree */ vmin = SDgetTreMin(st); - if (vmin <= FTINY) - return .0; + if ((vmin >= vmaxmin) | (vmin <= .01/M_PI)) + return .0; /* not worth bothering about */ SDsubtractTreVal(st, vmin); return M_PI * vmin; /* return hemispherical value */ } +/* Struct used in callback to find RGB extrema */ +typedef struct { + SDNode **stc; /* original Y, u' & v' trees */ + float rgb[3]; /* RGB value */ + SDNode *new_stu, *new_stv; /* replacement u' & v' trees */ +} SDextRGBs; + +/* Callback to find minimum RGB from Y value plus CIE (u',v') trees */ +static int +get_min_RGB(float yval, const double *cmin, double csiz, void *cptr) +{ + SDextRGBs *mp = (SDextRGBs *)cptr; + double cmax[SD_MAXDIM]; + float rgb[3]; + + if (mp->stc[tt_Y]->ndim == 3) { + if (cmin[0] + .5*csiz >= .5) + return 0; /* ignore dead half of isotropic */ + } else + cmax[3] = cmin[3] + csiz; + cmax[0] = cmin[0] + csiz; + cmax[1] = cmin[1] + csiz; + cmax[2] = cmin[2] + csiz; + /* average RGB color over voxel */ + SDyuv2rgb(yval, SDavgTreBox(mp->stc[tt_u], cmin, cmax), + SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb); + /* track smallest components */ + if (rgb[0] < mp->rgb[0]) mp->rgb[0] = rgb[0]; + if (rgb[1] < mp->rgb[1]) mp->rgb[1] = rgb[1]; + if (rgb[2] < mp->rgb[2]) mp->rgb[2] = rgb[2]; + return 0; +} + +/* Callback to build adjusted u' tree */ +static int +adjust_utree(float uprime, const double *cmin, double csiz, void *cptr) +{ + SDextRGBs *mp = (SDextRGBs *)cptr; + double cmax[SD_MAXDIM]; + double yval; + float rgb[3]; + C_COLOR clr; + + if (mp->stc[tt_Y]->ndim == 3) { + if (cmin[0] + .5*csiz >= .5) + return 0; /* ignore dead half of isotropic */ + } else + cmax[3] = cmin[3] + csiz; + cmax[0] = cmin[0] + csiz; + cmax[1] = cmin[1] + csiz; + cmax[2] = cmin[2] + csiz; + /* average RGB color over voxel */ + SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax), uprime, + SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb); + /* subtract minimum (& clamp) */ + if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval; + if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval; + if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval; + c_fromSharpRGB(rgb, &clr); /* compute new u' for adj. RGB */ + uprime = 4.*clr.cx/(-2.*clr.cx + 12.*clr.cy + 3.); + /* assign in new u' tree */ + mp->new_stu = SDsetVoxel(mp->new_stu, mp->stc[tt_Y]->ndim, + cmin, csiz, uprime); + return -(mp->new_stu == NULL); +} + +/* Callback to build adjusted v' tree */ +static int +adjust_vtree(float vprime, const double *cmin, double csiz, void *cptr) +{ + SDextRGBs *mp = (SDextRGBs *)cptr; + double cmax[SD_MAXDIM]; + double yval; + float rgb[3]; + C_COLOR clr; + + if (mp->stc[tt_Y]->ndim == 3) { + if (cmin[0] + .5*csiz >= .5) + return 0; /* ignore dead half of isotropic */ + } else + cmax[3] = cmin[3] + csiz; + cmax[0] = cmin[0] + csiz; + cmax[1] = cmin[1] + csiz; + cmax[2] = cmin[2] + csiz; + /* average RGB color over voxel */ + SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax), + SDavgTreBox(mp->stc[tt_u], cmin, cmax), + vprime, rgb); + /* subtract minimum (& clamp) */ + if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval; + if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval; + if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval; + c_fromSharpRGB(rgb, &clr); /* compute new v' for adj. RGB */ + vprime = 9.*clr.cy/(-2.*clr.cx + 12.*clr.cy + 3.); + /* assign in new v' tree */ + mp->new_stv = SDsetVoxel(mp->new_stv, mp->stc[tt_Y]->ndim, + cmin, csiz, vprime); + return -(mp->new_stv == NULL); +} + +/* Subtract minimum (diffuse) color and return luminance & CIE (x,y) */ +static double +subtract_min_RGB(C_COLOR *cs, SDNode *stc[]) +{ + SDextRGBs my_min; + double ymin; + + my_min.stc = stc; + my_min.rgb[0] = my_min.rgb[1] = my_min.rgb[2] = FHUGE; + my_min.new_stu = my_min.new_stv = NULL; + /* get minimum RGB value */ + SDtraverseTre(stc[tt_Y], NULL, 0, get_min_RGB, &my_min); + /* convert to C_COLOR */ + ymin = c_fromSharpRGB(my_min.rgb, cs); + if ((ymin >= .5*FHUGE) | (ymin <= .01/M_PI)) + return .0; /* close to zero or no tree */ + /* adjust u' & v' trees */ + SDtraverseTre(stc[tt_u], NULL, 0, adjust_utree, &my_min); + SDtraverseTre(stc[tt_v], NULL, 0, adjust_vtree, &my_min); + SDfreeTre(stc[tt_u]); SDfreeTre(stc[tt_v]); + stc[tt_u] = SDsimplifyTre(my_min.new_stu); + stc[tt_v] = SDsimplifyTre(my_min.new_stv); + /* subtract Y & return hemispherical */ + SDsubtractTreVal(stc[tt_Y], ymin); + + return M_PI * ymin; +} + /* Extract and separate diffuse portion of BSDF */ static void extract_diffuse(SDValue *dv, SDSpectralDF *df) { int n; + SDTre *sdt; if (df == NULL || df->ncomp <= 0) { dv->spec = c_dfcolor; dv->cieY = .0; return; } - dv->spec = df->comp[0].cspec[0]; - dv->cieY = subtract_min((*(SDTre *)df->comp[0].dist).st); - /* in case of multiple components */ - for (n = df->ncomp; --n; ) { - double ymin = subtract_min((*(SDTre *)df->comp[n].dist).st); - c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]); - dv->cieY += ymin; + sdt = (SDTre *)df->comp[0].dist; + /* subtract minimum color/grayscale */ + if (sdt->stc[tt_u] != NULL && sdt->stc[tt_v] != NULL) { + int i = 3*(tt_RGB_coef[1] < .001); + while (i--) { /* initialize on first call */ + float rgb[3]; + rgb[0] = rgb[1] = rgb[2] = .0f; rgb[i] = 1.f; + tt_RGB_coef[i] = c_fromSharpRGB(rgb, &tt_RGB_prim[i]); + } + memcpy(df->comp[0].cspec, tt_RGB_prim, sizeof(tt_RGB_prim)); + dv->cieY = subtract_min_RGB(&dv->spec, sdt->stc); + } else { + df->comp[0].cspec[0] = dv->spec = c_dfcolor; + dv->cieY = subtract_min_Y(sdt->stc[tt_Y]); } df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */ - /* make sure everything is set */ - c_ccvt(&dv->spec, C_CSXY+C_CSSPEC); + + c_ccvt(&dv->spec, C_CSXY); /* make sure (x,y) is set */ } /* Load a variable-resolution BSDF tree from an open XML file */ @@ -1047,24 +1395,34 @@ SDloadTre(SDData *sd, ezxml_t wtl) /* load BSDF components */ for (wld = ezxml_child(wtl, "WavelengthData"); wld != NULL; wld = wld->next) { - if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")), - "Visible")) - continue; /* just visible for now */ + const char *cnm = ezxml_txt(ezxml_child(wld,"Wavelength")); + int ct = -1; + if (!strcasecmp(cnm, "Visible")) + ct = tt_Y; + else if (!strcasecmp(cnm, "CIE-u")) + ct = tt_u; + else if (!strcasecmp(cnm, "CIE-v")) + ct = tt_v; + else + continue; for (wdb = ezxml_child(wld, "WavelengthDataBlock"); wdb != NULL; wdb = wdb->next) - if ((ec = load_bsdf_data(sd, wdb, rank)) != SDEnone) + if ((ec = load_bsdf_data(sd, wdb, ct, rank)) != SDEnone) return ec; } /* 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; } /* Variable resolution BSDF methods */ -SDFunc SDhandleTre = { +const SDFunc SDhandleTre = { &SDgetTreBSDF, &SDqueryTreProjSA, &SDgetTreCDist,