--- ray/src/cv/bsdfrep.c 2014/08/22 05:38:44 2.27 +++ ray/src/cv/bsdfrep.c 2021/12/15 01:38:50 2.37 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdfrep.c,v 2.27 2014/08/22 05:38:44 greg Exp $"; +static const char RCSid[] = "$Id: bsdfrep.c,v 2.37 2021/12/15 01:38:50 greg Exp $"; #endif /* * Support BSDF representation as radial basis functions. @@ -9,11 +9,11 @@ static const char RCSid[] = "$Id: bsdfrep.c,v 2.27 201 #define _USE_MATH_DEFINES #include -#include #include #include "rtio.h" #include "resolu.h" #include "bsdfrep.h" +#include "random.h" /* name and manufacturer if known */ char bsdf_name[256]; char bsdf_manuf[256]; @@ -29,12 +29,19 @@ int single_plane_incident = -1; int input_orient = 0; int output_orient = 0; + /* represented color space */ +RBColor rbf_colorimetry = RBCunknown; + +const char *RBCident[] = { + "CIE-Y", "CIE-XYZ", "Spectral", "Unknown" + }; + /* BSDF histogram */ unsigned long bsdf_hist[HISTLEN]; /* BSDF value for boundary regions */ double bsdf_min = 0; -double bsdf_spec_peak = 0; +double bsdf_spec_val = 0; double bsdf_spec_rad = 0; /* processed incident DSF measurements */ @@ -46,18 +53,13 @@ MIGRATION *mig_list = NULL; /* current input direction */ double theta_in_deg, phi_in_deg; + /* header line sharing callback */ +int (*sir_headshare)(char *s) = NULL; + /* Register new input direction */ int new_input_direction(double new_theta, double new_phi) { - if (!input_orient) /* check input orientation */ - input_orient = 1 - 2*(new_theta > 90.); - else if (input_orient > 0 ^ new_theta < 90.) { - fprintf(stderr, - "%s: Cannot handle input angles on both sides of surface\n", - progname); - return(0); - } /* normalize angle ranges */ while (new_theta < -180.) new_theta += 360.; @@ -67,12 +69,21 @@ new_input_direction(double new_theta, double new_phi) new_theta = -new_theta; new_phi += 180.; } - if ((theta_in_deg = new_theta) < 1.0) - return(1); /* don't rely on phi near normal */ while (new_phi < 0) new_phi += 360.; while (new_phi >= 360.) new_phi -= 360.; + /* check input orientation */ + if (!input_orient) + input_orient = 1 - 2*(new_theta > 90.); + else if (input_orient > 0 ^ new_theta < 90.) { + fprintf(stderr, + "%s: Cannot handle input angles on both sides of surface\n", + progname); + return(0); + } + if ((theta_in_deg = new_theta) < 1.0) + return(1); /* don't rely on phi near normal */ if (single_plane_incident > 0) /* check input coverage */ single_plane_incident = (round(new_phi) == round(phi_in_deg)); else if (single_plane_incident < 0) @@ -93,7 +104,10 @@ new_input_direction(double new_theta, double new_phi) int use_symmetry(FVECT vec) { - const double phi = get_phi360(vec); + double phi = get_phi360(vec); + /* because of -0. issue */ + while (phi >= 360.) phi -= 360.; + while (phi < 0.) phi += 360.; switch (inp_coverage) { case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4: @@ -214,10 +228,10 @@ rotate_rbf(RBFNODE *rbf, const FVECT invec) void ovec_from_pos(FVECT vec, int xpos, int ypos) { - double uv[2]; + RREAL uv[2]; double r2; - SDsquare2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res); + square2disk(uv, (xpos+.5)/grid_res, (ypos+.5)/grid_res); /* uniform hemispherical projection */ r2 = uv[0]*uv[0] + uv[1]*uv[1]; vec[0] = vec[1] = sqrt(2. - r2); @@ -230,10 +244,10 @@ ovec_from_pos(FVECT vec, int xpos, int ypos) void pos_from_vec(int pos[2], const FVECT vec) { - double sq[2]; /* uniform hemispherical projection */ + RREAL sq[2]; /* uniform hemispherical projection */ double norm = 1./sqrt(1. + fabs(vec[2])); - SDdisk2square(sq, vec[0]*norm, vec[1]*norm); + disk2square(sq, vec[0]*norm, vec[1]*norm); pos[0] = (int)(sq[0]*grid_res); pos[1] = (int)(sq[1]*grid_res); @@ -259,23 +273,28 @@ rbf_volume(const RBFVAL *rbfp) return(integ); } -/* Evaluate BSDF at the given normalized outgoing direction */ -double -eval_rbfrep(const RBFNODE *rp, const FVECT outvec) +/* Evaluate BSDF at the given normalized outgoing direction in color */ +SDError +eval_rbfcol(SDValue *sv, const RBFNODE *rp, const FVECT outvec) { const double rfact2 = (38./M_PI/M_PI)*(grid_res*grid_res); int pos[2]; double res = 0; + double usum = 0, vsum = 0; const RBFVAL *rbfp; FVECT odir; double rad2; int n; + /* assign default value */ + sv->spec = c_dfcolor; + sv->cieY = bsdf_min; /* check for wrong side */ - if (outvec[2] > 0 ^ output_orient > 0) - return(.0); - /* use minimum if no information avail. */ - if (rp == NULL) - return(bsdf_min); + if (outvec[2] > 0 ^ output_orient > 0) { + strcpy(SDerrorDetail, "Wrong-side scattering query"); + return(SDEargument); + } + if (rp == NULL) /* return minimum if no information avail. */ + return(SDEnone); /* optimization for fast lobe culling */ pos_from_vec(pos, outvec); /* sum radial basis function */ @@ -283,19 +302,42 @@ eval_rbfrep(const RBFNODE *rp, const FVECT outvec) for (n = rp->nrbf; n--; rbfp++) { int d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) + (pos[1]-rbfp->gy)*(pos[1]-rbfp->gy); + double val; rad2 = R2ANG(rbfp->crad); rad2 *= rad2; if (d2 > rad2*rfact2) continue; ovec_from_pos(odir, rbfp->gx, rbfp->gy); - res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2); + val = rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2); + if (rbf_colorimetry == RBCtristimulus) { + usum += val * (rbfp->chroma & 0xff); + vsum += val * (rbfp->chroma>>8 & 0xff); + } + res += val; } - res /= COSF(outvec[2]); - if (res < bsdf_min) /* never return less than bsdf_min */ - return(bsdf_min); - return(res); + sv->cieY = res / COSF(outvec[2]); + if (sv->cieY < bsdf_min) { /* never return less than bsdf_min */ + sv->cieY = bsdf_min; + } else if (rbf_colorimetry == RBCtristimulus) { + C_CHROMA cres = (int)(usum/res + frandom()); + cres |= (int)(vsum/res + frandom()) << 8; + c_decodeChroma(&sv->spec, cres); + } + return(SDEnone); } +/* Evaluate BSDF at the given normalized outgoing direction in Y */ +double +eval_rbfrep(const RBFNODE *rp, const FVECT outvec) +{ + SDValue sv; + + if (eval_rbfcol(&sv, rp, outvec) == SDEnone) + return(sv.cieY); + + return(0.0); +} + /* Insert a new directional scattering function in our global list */ int insert_dsf(RBFNODE *newrbf) @@ -406,7 +448,7 @@ def_rbf_spec(const FVECT invec) if (input_orient > 0 ^ invec[2] > 0) /* wrong side? */ return(NULL); - if ((bsdf_spec_peak <= bsdf_min) | (bsdf_spec_rad <= 0)) + if ((bsdf_spec_val <= bsdf_min) | (bsdf_spec_rad <= 0)) return(NULL); /* nothing set */ rbf = (RBFNODE *)malloc(sizeof(RBFNODE)); if (rbf == NULL) @@ -420,7 +462,8 @@ def_rbf_spec(const FVECT invec) rbf->ejl = NULL; VCOPY(rbf->invec, invec); rbf->nrbf = 1; - rbf->rbfa[0].peak = bsdf_spec_peak * output_orient*ovec[2]; + rbf->rbfa[0].peak = bsdf_spec_val * COSF(ovec[2]); + rbf->rbfa[0].chroma = c_dfchroma; rbf->rbfa[0].crad = ANG2R(bsdf_spec_rad); rbf->rbfa[0].gx = pos[0]; rbf->rbfa[0].gy = pos[1]; @@ -438,7 +481,7 @@ e_advect_rbf(const MIGRATION *mig, const FVECT invec, double t, full_dist; /* get relative position */ t = Acos(DOT(invec, mig->rbfv[0]->invec)); - if (t < M_PI/grid_res) { /* near first DSF */ + if (t <= .001) { /* near first DSF */ n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1); rbf = (RBFNODE *)malloc(n); if (rbf == NULL) @@ -448,7 +491,7 @@ e_advect_rbf(const MIGRATION *mig, const FVECT invec, return(rbf); } full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec)); - if (t > full_dist-M_PI/grid_res) { /* near second DSF */ + if (t >= full_dist-.001) { /* near second DSF */ n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1); rbf = (RBFNODE *)malloc(n); if (rbf == NULL) @@ -484,9 +527,11 @@ tryagain: const RBFVAL *rbf0i = &mig->rbfv[0]->rbfa[i]; const float peak0 = rbf0i->peak; const double rad0 = R2ANG(rbf0i->crad); + C_COLOR cc0; FVECT v0; float mv; ovec_from_pos(v0, rbf0i->gx, rbf0i->gy); + c_decodeChroma(&cc0, rbf0i->chroma); for (j = 0; j < mtx_ncols(mig); j++) if ((mv = mtx_coef(mig,i,j)) > cthresh) { const RBFVAL *rbf1j = &mig->rbfv[1]->rbfa[j]; @@ -497,6 +542,13 @@ tryagain: rad2 = rad0*rad0*(1.-t) + rad2*rad2*t; rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal * rad0*rad0/rad2; + if (rbf_colorimetry == RBCtristimulus) { + C_COLOR cres; + c_decodeChroma(&cres, rbf1j->chroma); + c_cmix(&cres, 1.-t, &cc0, t, &cres); + rbf->rbfa[n].chroma = c_encodeChroma(&cres); + } else + rbf->rbfa[n].chroma = c_dfchroma; rbf->rbfa[n].crad = ANG2R(sqrt(rad2)); ovec_from_pos(v, rbf1j->gx, rbf1j->gy); geodesic(v, v0, v, t, GEOD_REL); @@ -533,9 +585,11 @@ clear_bsdf_rep(void) inp_coverage = 0; single_plane_incident = -1; input_orient = output_orient = 0; + rbf_colorimetry = RBCunknown; grid_res = GRIDRES; + memset(bsdf_hist, 0, sizeof(bsdf_hist)); bsdf_min = 0; - bsdf_spec_peak = 0; + bsdf_spec_val = 0; bsdf_spec_rad = 0; } @@ -553,12 +607,14 @@ save_bsdf_rep(FILE *ofp) fprintf(ofp, "MANUFACT=%s\n", bsdf_manuf); fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage); fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient); + fprintf(ofp, "COLORIMETRY=%s\n", RBCident[rbf_colorimetry]); fprintf(ofp, "GRIDRES=%d\n", grid_res); fprintf(ofp, "BSDFMIN=%g\n", bsdf_min); - if ((bsdf_spec_peak > bsdf_min) & (bsdf_spec_rad > 0)) - fprintf(ofp, "BSDFSPEC= %f %f\n", bsdf_spec_peak, bsdf_spec_rad); + if ((bsdf_spec_val > bsdf_min) & (bsdf_spec_rad > 0)) + fprintf(ofp, "BSDFSPEC= %f %f\n", bsdf_spec_val, bsdf_spec_rad); fputformat(BSDFREP_FMT, ofp); fputc('\n', ofp); + putint(BSDFREP_MAGIC, 2, ofp); /* write each DSF */ for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) { putint(rbf->ord, 4, ofp); @@ -569,9 +625,10 @@ save_bsdf_rep(FILE *ofp) putint(rbf->nrbf, 4, ofp); for (i = 0; i < rbf->nrbf; i++) { putflt(rbf->rbfa[i].peak, ofp); + putint(rbf->rbfa[i].chroma, 2, ofp); putint(rbf->rbfa[i].crad, 2, ofp); - putint(rbf->rbfa[i].gx, 1, ofp); - putint(rbf->rbfa[i].gy, 1, ofp); + putint(rbf->rbfa[i].gx, 2, ofp); + putint(rbf->rbfa[i].gy, 2, ofp); } } putint(-1, 4, ofp); /* terminator */ @@ -607,39 +664,57 @@ save_bsdf_rep(FILE *ofp) static int headline(char *s, void *p) { - char fmt[32]; + char fmt[MAXFMTLEN]; + int i; + if (isheadid(s)) + return(0); if (!strncmp(s, "NAME=", 5)) { strcpy(bsdf_name, s+5); bsdf_name[strlen(bsdf_name)-1] = '\0'; + return(1); } if (!strncmp(s, "MANUFACT=", 9)) { strcpy(bsdf_manuf, s+9); bsdf_manuf[strlen(bsdf_manuf)-1] = '\0'; + return(1); } if (!strncmp(s, "SYMMETRY=", 9)) { inp_coverage = atoi(s+9); single_plane_incident = !inp_coverage; - return(0); + return(1); } if (!strncmp(s, "IO_SIDES=", 9)) { sscanf(s+9, "%d %d", &input_orient, &output_orient); - return(0); + return(1); } + if (!strncmp(s, "COLORIMETRY=", 12)) { + fmt[0] = '\0'; + sscanf(s+12, "%s", fmt); + for (i = RBCunknown; i >= 0; i--) + if (!strcmp(fmt, RBCident[i])) + break; + if (i < 0) + return(-1); + rbf_colorimetry = i; + return(1); + } if (!strncmp(s, "GRIDRES=", 8)) { sscanf(s+8, "%d", &grid_res); - return(0); + return(1); } if (!strncmp(s, "BSDFMIN=", 8)) { sscanf(s+8, "%lf", &bsdf_min); - return(0); + return(1); } if (!strncmp(s, "BSDFSPEC=", 9)) { - sscanf(s+9, "%lf %lf", &bsdf_spec_peak, &bsdf_spec_rad); - return(0); + sscanf(s+9, "%lf %lf", &bsdf_spec_val, &bsdf_spec_rad); + return(1); } - if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT)) - return(-1); + if (formatval(fmt, s)) + return (strcmp(fmt, BSDFREP_FMT) ? -1 : 0); + if (sir_headshare != NULL) + return ((*sir_headshare)(s)); return(0); } @@ -656,11 +731,16 @@ load_bsdf_rep(FILE *ifp) return(0); if (getheader(ifp, headline, NULL) < 0 || (single_plane_incident < 0) | !input_orient | !output_orient | - (grid_res < 16) | (grid_res > 256)) { + (grid_res < 16) | (grid_res > 0xffff)) { fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n", progname); return(0); } + if (getint(2, ifp) != BSDFREP_MAGIC) { + fprintf(stderr, "%s: bad magic number for BSDF interpolant\n", + progname); + return(0); + } memset(&rbfh, 0, sizeof(rbfh)); /* read each DSF */ while ((rbfh.ord = getint(4, ifp)) >= 0) { RBFNODE *newrbf; @@ -681,9 +761,10 @@ load_bsdf_rep(FILE *ifp) *newrbf = rbfh; for (i = 0; i < rbfh.nrbf; i++) { newrbf->rbfa[i].peak = getflt(ifp); + newrbf->rbfa[i].chroma = getint(2, ifp) & 0xffff; newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff; - newrbf->rbfa[i].gx = getint(1, ifp) & 0xff; - newrbf->rbfa[i].gy = getint(1, ifp) & 0xff; + newrbf->rbfa[i].gx = getint(2, ifp) & 0xffff; + newrbf->rbfa[i].gy = getint(2, ifp) & 0xffff; } if (feof(ifp)) goto badEOF;