--- ray/src/cv/bsdf2klems.c 2016/02/03 01:57:06 2.20 +++ ray/src/cv/bsdf2klems.c 2020/10/26 21:12:20 2.29 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf2klems.c,v 2.20 2016/02/03 01:57:06 greg Exp $"; +static const char RCSid[] = "$Id: bsdf2klems.c,v 2.29 2020/10/26 21:12:20 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.20 */ #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) { @@ -202,15 +201,14 @@ eval_bsdf(const char *fname) 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; + xsum += sdv.cieY * sdv.spec.cx; + ysum += sdv.cieY * sdv.spec.cy; } } fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps); if (rbf_colorimetry == RBCtristimulus) { - fprintf(cfp[CIE_X], "\t%3e\n", xsum*sum/(npsamps*ysum)); - fprintf(cfp[CIE_Z], "\t%3e\n", + fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum)); + fprintf(cfp[CIE_Z], "\t%.3e\n", (sum - xsum - ysum)*sum/(npsamps*ysum)); } } @@ -252,15 +250,14 @@ eval_bsdf(const char *fname) 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; + xsum += sdv.cieY * sdv.spec.cx; + ysum += sdv.cieY * sdv.spec.cy; } } fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps); if (rbf_colorimetry == RBCtristimulus) { - fprintf(cfp[CIE_X], "\t%3e\n", xsum*sum/(npsamps*ysum)); - fprintf(cfp[CIE_Z], "\t%3e\n", + fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum)); + fprintf(cfp[CIE_Z], "\t%.3e\n", (sum - xsum - ysum)*sum/(npsamps*ysum)); } } @@ -301,15 +298,14 @@ eval_bsdf(const char *fname) 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; + xsum += sdv.cieY * sdv.spec.cx; + ysum += sdv.cieY * sdv.spec.cy; } } fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps); if (rbf_colorimetry == RBCtristimulus) { - fprintf(cfp[CIE_X], "\t%3e\n", xsum*sum/(npsamps*ysum)); - fprintf(cfp[CIE_Z], "\t%3e\n", + fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum)); + fprintf(cfp[CIE_Z], "\t%.3e\n", (sum - xsum - ysum)*sum/(npsamps*ysum)); } } @@ -332,12 +328,13 @@ eval_bsdf(const char *fname) if ((bsd.tb != NULL) | (bsd.tf != NULL)) { input_orient = -1; output_orient = 1; cfp[CIE_Y] = open_component_file(CIE_Y); - if (bsd.tb != NULL && bsd.tb->comp[0].cspec[2].flags) { - rbf_colorimetry = RBCtristimulus; + if (bsd.tb != NULL) + rbf_colorimetry = bsd.tb->comp[0].cspec[2].flags + ? RBCtristimulus : RBCphotopic ; + if (rbf_colorimetry == RBCtristimulus) { cfp[CIE_X] = open_component_file(CIE_X); cfp[CIE_Z] = open_component_file(CIE_Z); - } else - rbf_colorimetry = RBCphotopic; + } for (j = 0; j < abp->nangles; j++) { for (i = 0; i < abp->nangles; i++) { sum = 0; /* average over patches */ @@ -350,15 +347,14 @@ eval_bsdf(const char *fname) 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; + xsum += sdv.cieY * sdv.spec.cx; + ysum += sdv.cieY * sdv.spec.cy; } } fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps); if (rbf_colorimetry == RBCtristimulus) { - fprintf(cfp[CIE_X], "\t%3e\n", xsum*sum/(npsamps*ysum)); - fprintf(cfp[CIE_Z], "\t%3e\n", + fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum)); + fprintf(cfp[CIE_Z], "\t%.3e\n", (sum - xsum - ysum)*sum/(npsamps*ysum)); } } @@ -434,92 +430,97 @@ 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); - xsum += sdv.cieY*sdv.spec.cx; - ysum += sdv.cieY*sdv.spec.cy; - } + if (rbf_colorimetry == RBCtristimulus) { + xsum += sdv.cieY * sdv.spec.cx; + ysum += sdv.cieY * sdv.spec.cy; + } } - n = j*abp->nangles + i; - bsdfarr[n] = sum / (double)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) @@ -575,6 +576,39 @@ wrap_up(void) } #endif +#define HEAD_BUFLEN 8192 +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; + + 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[]) @@ -606,7 +640,15 @@ main(int argc, char *argv[]) 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 @@ -691,9 +733,13 @@ main(int argc, char *argv[]) progname, argv[i]); return(1); } + sprintf(pbuf, "%s:\n", argv[i]); + record2header(pbuf); + sir_headshare = &record2header; if (!load_bsdf_rep(fpin)) return(1); fclose(fpin); + done_header(); sprintf(pbuf, "Interpolating component '%s'", argv[i]); prog_start(pbuf); eval_rbf();