--- ray/src/cv/bsdf2klems.c 2016/02/04 00:45:47 2.21 +++ ray/src/cv/bsdf2klems.c 2023/07/26 15:50:03 2.35 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf2klems.c,v 2.21 2016/02/04 00:45:47 greg Exp $"; +static const char RCSid[] = "$Id: bsdf2klems.c,v 2.35 2023/07/26 15:50:03 greg Exp $"; #endif /* * Load measured BSDF interpolant and write out as XML file with Klems matrix. @@ -8,10 +8,9 @@ static const char RCSid[] = "$Id: bsdf2klems.c,v 2.21 */ #define _USE_MATH_DEFINES -#include #include -#include #include +#include #include "random.h" #include "platform.h" #include "paths.h" @@ -31,7 +30,7 @@ static const char klems_half[] = "LBNL/Klems Half"; static const char klems_quarter[] = "LBNL/Klems Quarter"; static const char *kbasis = klems_full; /* number of BSDF samples per patch */ -static int npsamps = 256; +static int npsamps = 1024; /* limit on number of RBF lobes */ static int lobe_lim = 15000; /* progress bar length */ @@ -41,7 +40,7 @@ static int do_prog = 79; static char *wrapBSDF[MAXCARG] = {"wrapBSDF", "-W", "-UU"}; static int wbsdfac = 3; -/* Add argument to wrapBSDF, allocating space if isstatic */ +/* Add argument to wrapBSDF, allocating space if !isstatic */ static void add_wbsdf(const char *arg, int isstatic) { @@ -197,12 +196,11 @@ eval_bsdf(const char *fname) for (n = npsamps; n-- > 0; ) { fo_getvec(vout, j+(n+frandom())/npsamps, abp); fi_getvec(vin, i+urand(n), abp); - ec = SDevalBSDF(&sdv, vout, vin, &bsd); + ec = SDevalBSDF(&sdv, vin, vout, &bsd); if (ec != SDEnone) goto err; sum += sdv.cieY; if (rbf_colorimetry == RBCtristimulus) { - c_ccvt(&sdv.spec, C_CSXY); xsum += sdv.cieY * sdv.spec.cx; ysum += sdv.cieY * sdv.spec.cy; } @@ -247,12 +245,11 @@ eval_bsdf(const char *fname) for (n = npsamps; n-- > 0; ) { bo_getvec(vout, j+(n+frandom())/npsamps, abp); bi_getvec(vin, i+urand(n), abp); - ec = SDevalBSDF(&sdv, vout, vin, &bsd); + ec = SDevalBSDF(&sdv, vin, vout, &bsd); if (ec != SDEnone) goto err; sum += sdv.cieY; if (rbf_colorimetry == RBCtristimulus) { - c_ccvt(&sdv.spec, C_CSXY); xsum += sdv.cieY * sdv.spec.cx; ysum += sdv.cieY * sdv.spec.cy; } @@ -280,7 +277,7 @@ eval_bsdf(const char *fname) } } /* front transmission */ - if (bsd.tf != NULL || bsd.tLamb.cieY > .002) { + if (bsd.tf != NULL || bsd.tLambFront.cieY > .002) { input_orient = 1; output_orient = -1; cfp[CIE_Y] = open_component_file(CIE_Y); if (bsd.tf != NULL && bsd.tf->comp[0].cspec[2].flags) { @@ -296,12 +293,11 @@ eval_bsdf(const char *fname) for (n = npsamps; n-- > 0; ) { bo_getvec(vout, j+(n+frandom())/npsamps, abp); fi_getvec(vin, i+urand(n), abp); - ec = SDevalBSDF(&sdv, vout, vin, &bsd); + ec = SDevalBSDF(&sdv, vin, vout, &bsd); if (ec != SDEnone) goto err; sum += sdv.cieY; if (rbf_colorimetry == RBCtristimulus) { - c_ccvt(&sdv.spec, C_CSXY); xsum += sdv.cieY * sdv.spec.cx; ysum += sdv.cieY * sdv.spec.cy; } @@ -346,12 +342,11 @@ eval_bsdf(const char *fname) for (n = npsamps; n-- > 0; ) { fo_getvec(vout, j+(n+frandom())/npsamps, abp); bi_getvec(vin, i+urand(n), abp); - ec = SDevalBSDF(&sdv, vout, vin, &bsd); + ec = SDevalBSDF(&sdv, vin, vout, &bsd); if (ec != SDEnone) goto err; sum += sdv.cieY; if (rbf_colorimetry == RBCtristimulus) { - c_ccvt(&sdv.spec, C_CSXY); xsum += sdv.cieY * sdv.spec.cx; ysum += sdv.cieY * sdv.spec.cy; } @@ -435,97 +430,102 @@ eval_function(char *funame) static void eval_rbf(void) { - ANGLE_BASIS *abp = get_basis(kbasis); - float (*XZarr)[2] = NULL; - float bsdfarr[MAXPATCHES*MAXPATCHES]; - FILE *cfp[3]; - FVECT vin, vout; - double sum, xsum, ysum; - int i, j, n; - /* sanity check */ - if (abp->nangles > MAXPATCHES) { - fprintf(stderr, "%s: too many patches!\n", progname); - exit(1); - } - if (rbf_colorimetry == RBCtristimulus) - XZarr = (float (*)[2])malloc(sizeof(float)*2*abp->nangles*abp->nangles); - for (i = 0; i < abp->nangles; i++) { - RBFNODE *rbf; - if (input_orient > 0) /* use incident patch center */ - fi_getvec(vin, i+.5*(i>0), abp); - else - bi_getvec(vin, i+.5*(i>0), abp); + ANGLE_BASIS *abp = get_basis(kbasis); + float (*XZarr)[2] = NULL; + float bsdfarr[MAXPATCHES*MAXPATCHES]; + FILE *cfp[3]; + FVECT vin, vout; + double sum, xsum, ysum, normf; + int i, j, ni, no, nisamps, nosamps; + /* sanity check */ + if (abp->nangles > MAXPATCHES) { + fprintf(stderr, "%s: too many patches!\n", progname); + exit(1); + } + memset(bsdfarr, 0, sizeof(bsdfarr)); + if (rbf_colorimetry == RBCtristimulus) + XZarr = (float (*)[2])calloc(abp->nangles*abp->nangles, 2*sizeof(float)); + nosamps = (int)(pow((double)npsamps, 0.67) + .5); + nisamps = (npsamps + (nosamps>>1)) / nosamps; + normf = 1./(double)(nisamps*nosamps); + for (i = 0; i < abp->nangles; i++) { + for (ni = nisamps; ni--; ) { /* sample over incident patch */ + RBFNODE *rbf; + if (input_orient > 0) /* vary incident patch loc. */ + fi_getvec(vin, i+urand(ni), abp); + else + bi_getvec(vin, i+urand(ni), abp); - rbf = advect_rbf(vin, lobe_lim); /* compute radial basis func */ + rbf = advect_rbf(vin, lobe_lim); /* compute radial basis func */ - for (j = 0; j < abp->nangles; j++) { - sum = 0; /* sample over exiting patch */ + for (j = 0; j < abp->nangles; j++) { + sum = 0; /* sample over exiting patch */ xsum = ysum = 0; - for (n = npsamps; n--; ) { + for (no = nosamps; no--; ) { SDValue sdv; if (output_orient > 0) - fo_getvec(vout, j+(n+frandom())/npsamps, abp); + fo_getvec(vout, j+(no+frandom())/nosamps, abp); else - bo_getvec(vout, j+(n+frandom())/npsamps, abp); + bo_getvec(vout, j+(no+frandom())/nosamps, abp); eval_rbfcol(&sdv, rbf, vout); sum += sdv.cieY; - if (XZarr != NULL) { - c_ccvt(&sdv.spec, C_CSXY); + if (rbf_colorimetry == RBCtristimulus) { xsum += sdv.cieY * sdv.spec.cx; ysum += sdv.cieY * sdv.spec.cy; - } + } } - n = j*abp->nangles + i; - bsdfarr[n] = sum / npsamps; - if (XZarr != NULL) { - XZarr[n][0] = xsum*sum/(npsamps*ysum); - XZarr[n][1] = (sum - xsum - ysum)*sum/(npsamps*ysum); + no = j*abp->nangles + i; + bsdfarr[no] += sum * normf; + if (rbf_colorimetry == RBCtristimulus) { + XZarr[no][0] += xsum*sum*normf/ysum; + XZarr[no][1] += (sum - xsum - ysum)*sum*normf/ysum; } } - if (rbf != NULL) + if (rbf != NULL) free(rbf); - prog_show((i+1.)/abp->nangles); } - /* write out our matrix */ - cfp[CIE_Y] = open_component_file(CIE_Y); - n = 0; - for (j = 0; j < abp->nangles; j++) { - for (i = 0; i < abp->nangles; i++, n++) - fprintf(cfp[CIE_Y], "\t%.3e\n", bsdfarr[n]); - fputc('\n', cfp[CIE_Y]); - } - prog_done(); - if (fclose(cfp[CIE_Y])) { - fprintf(stderr, "%s: error writing Y output\n", progname); - exit(1); - } - if (XZarr == NULL) /* no color? */ - return; - cfp[CIE_X] = open_component_file(CIE_X); - cfp[CIE_Z] = open_component_file(CIE_Z); - n = 0; - for (j = 0; j < abp->nangles; j++) { - for (i = 0; i < abp->nangles; i++, n++) { - fprintf(cfp[CIE_X], "\t%.3e\n", XZarr[n][0]); - fprintf(cfp[CIE_Z], "\t%.3e\n", XZarr[n][1]); - } - fputc('\n', cfp[CIE_X]); - fputc('\n', cfp[CIE_Z]); - } - free(XZarr); - if (fclose(cfp[CIE_X]) || fclose(cfp[CIE_Z])) { - fprintf(stderr, "%s: error writing X/Z output\n", progname); - exit(1); - } + prog_show((i+1.)/abp->nangles); + } + /* write out our matrix */ + cfp[CIE_Y] = open_component_file(CIE_Y); + no = 0; + for (j = 0; j < abp->nangles; j++) { + for (i = 0; i < abp->nangles; i++, no++) + fprintf(cfp[CIE_Y], "\t%.3e\n", bsdfarr[no]); + fputc('\n', cfp[CIE_Y]); + } + prog_done(); + if (fclose(cfp[CIE_Y])) { + fprintf(stderr, "%s: error writing Y output\n", progname); + exit(1); + } + if (XZarr == NULL) /* no color? */ + return; + cfp[CIE_X] = open_component_file(CIE_X); + cfp[CIE_Z] = open_component_file(CIE_Z); + no = 0; + for (j = 0; j < abp->nangles; j++) { + for (i = 0; i < abp->nangles; i++, no++) { + fprintf(cfp[CIE_X], "\t%.3e\n", XZarr[no][0]); + fprintf(cfp[CIE_Z], "\t%.3e\n", XZarr[no][1]); + } + fputc('\n', cfp[CIE_X]); + fputc('\n', cfp[CIE_Z]); + } + free(XZarr); + if (fclose(cfp[CIE_X]) || fclose(cfp[CIE_Z])) { + fprintf(stderr, "%s: error writing X/Z output\n", progname); + exit(1); + } } -#ifdef _WIN32 +#if defined(_WIN32) || defined(_WIN64) /* Execute wrapBSDF command (may never return) */ static int wrap_up(void) { - char cmd[8192]; + char cmd[32700]; if (bsdf_manuf[0]) { add_wbsdf("-f", 1); @@ -576,12 +576,49 @@ wrap_up(void) } #endif +#define HEAD_BUFLEN 10240 +static char head_buf[HEAD_BUFLEN]; +static int cur_headlen = 0; + +/* Record header line as comment associated with this SIR input */ +static int +record2header(char *s) +{ + int len = strlen(s); + + if (cur_headlen+len >= HEAD_BUFLEN-6) + return(0); + /* includes EOL */ + strcpy(head_buf+cur_headlen, s); + cur_headlen += len; + +#if defined(_WIN32) || defined(_WIN64) + if (head_buf[cur_headlen-1] == '\n') + head_buf[cur_headlen-1] = '\t'; +#endif + return(1); +} + +/* Finish off header for this file */ +static void +done_header(void) +{ + while (cur_headlen > 0 && isspace(head_buf[cur_headlen-1])) + --cur_headlen; + head_buf[cur_headlen] = '\0'; + if (!cur_headlen) + return; + add_wbsdf("-C", 1); + add_wbsdf(head_buf, 0); + head_buf[cur_headlen=0] = '\0'; +} + /* Read in BSDF and interpolate as Klems matrix representation */ int main(int argc, char *argv[]) { int dofwd = 0, dobwd = 1; - char buf[2048]; + char buf[1024]; char *cp; int i, na; @@ -602,12 +639,20 @@ main(int argc, char *argv[]) single_plane_incident = 0; break; case 'f': - if (!argv[i][2]) { + if ((argv[i][0] == '-') & !argv[i][2]) { if (strchr(argv[++i], '=') != NULL) { add_wbsdf("-f", 1); add_wbsdf(argv[i], 1); } else { - fcompile(argv[i]); + char *fpath = getpath(argv[i], + getrlibpath(), 0); + if (fpath == NULL) { + fprintf(stderr, + "%s: cannot find file '%s'\n", + argv[0], argv[i]); + return(1); + } + fcompile(fpath); single_plane_incident = 0; } } else @@ -685,25 +730,31 @@ main(int argc, char *argv[]) if (i < argc) { /* open input files if given */ int nbsdf = 0; for ( ; i < argc; i++) { /* interpolate each component */ - char pbuf[256]; FILE *fpin = fopen(argv[i], "rb"); if (fpin == NULL) { fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n", progname, argv[i]); return(1); } + sprintf(buf, "%s:\n", argv[i]); + record2header(buf); + sir_headshare = &record2header; if (!load_bsdf_rep(fpin)) return(1); fclose(fpin); - sprintf(pbuf, "Interpolating component '%s'", argv[i]); - prog_start(pbuf); + done_header(); + sprintf(buf, "Interpolating component '%s'", argv[i]); + prog_start(buf); eval_rbf(); } return(wrap_up()); } SET_FILE_BINARY(stdin); /* load from stdin */ + record2header(":\n"); + sir_headshare = &record2header; if (!load_bsdf_rep(stdin)) return(1); + done_header(); prog_start("Interpolating from standard input"); eval_rbf(); /* resample dist. */ return(wrap_up());