--- ray/src/cv/bsdf2klems.c 2018/05/04 23:56:49 2.25 +++ 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.25 2018/05/04 23:56:49 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.25 */ #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 */ @@ -431,88 +430,94 @@ 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 (rbf_colorimetry == RBCtristimulus) { xsum += sdv.cieY * sdv.spec.cx; ysum += sdv.cieY * sdv.spec.cy; - } + } } - n = j*abp->nangles + i; - bsdfarr[n] = sum / npsamps; + no = j*abp->nangles + i; + bsdfarr[no] += sum * normf; if (rbf_colorimetry == RBCtristimulus) { - XZarr[n][0] = xsum*sum/(npsamps*ysum); - XZarr[n][1] = (sum - xsum - ysum)*sum/(npsamps*ysum); + 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); + } } #if defined(_WIN32) || defined(_WIN64) @@ -571,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[]) @@ -695,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();