--- ray/src/cv/bsdf2ttree.c 2016/03/06 01:13:17 2.36 +++ ray/src/cv/bsdf2ttree.c 2020/05/08 22:34:21 2.49 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf2ttree.c,v 2.36 2016/03/06 01:13:17 schorsch Exp $"; +static const char RCSid[] = "$Id: bsdf2ttree.c,v 2.49 2020/05/08 22:34:21 greg Exp $"; #endif /* * Load measured BSDF interpolant and write out as XML file with tensor tree. @@ -8,7 +8,6 @@ static const char RCSid[] = "$Id: bsdf2ttree.c,v 2.36 */ #define _USE_MATH_DEFINES -#include #include #include #include "random.h" @@ -19,15 +18,17 @@ static const char RCSid[] = "$Id: bsdf2ttree.c,v 2.36 #include "bsdfrep.h" /* global argv[0] */ char *progname; + /* reciprocity averaging option */ +static const char *recip = " -a"; /* percentage to cull (<0 to turn off) */ static double pctcull = 90.; /* sampling order */ static int samp_order = 6; /* super-sampling threshold */ -const double ssamp_thresh = 0.35; +static double ssamp_thresh = 0.35; /* number of super-samples */ #ifndef NSSAMP -#define NSSAMP 100 +#define NSSAMP 256 #endif /* limit on number of RBF lobes */ static int lobe_lim = 15000; @@ -38,7 +39,7 @@ static int do_prog = 79; static char *wrapBSDF[MAXCARG] = {"wrapBSDF", "-U"}; static int wbsdfac = 2; -/* 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) { @@ -127,6 +128,9 @@ static void eval_isotropic(char *funame) { const int sqres = 1<= 0) { - sprintf(cmd, "rttree_reduce -a -h -ff -r 3 -t %f -g %d > %s", - pctcull, samp_order, create_component_file(0)); + sprintf(cmd, "rttree_reduce%s -h -ff -r 3 -t %f -g %d > %s", + recip, pctcull, samp_order, create_component_file(0)); ofp = popen(cmd, "w"); if (ofp == NULL) { fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n", @@ -149,11 +153,11 @@ eval_isotropic(char *funame) #endif if (rbf_colorimetry == RBCtristimulus) { double uvcull = 100. - (100.-pctcull)*.25; - sprintf(cmd, "rttree_reduce -a -h -ff -r 3 -t %f -g %d > %s", - uvcull, samp_order, create_component_file(1)); + sprintf(cmd, "rttree_reduce%s -h -ff -r 3 -t %f -g %d > %s", + recip, uvcull, samp_order, create_component_file(1)); uvfp[0] = popen(cmd, "w"); - sprintf(cmd, "rttree_reduce -a -h -ff -r 3 -t %f -g %d > %s", - uvcull, samp_order, create_component_file(2)); + sprintf(cmd, "rttree_reduce%s -h -ff -r 3 -t %f -g %d > %s", + recip, uvcull, samp_order, create_component_file(2)); uvfp[1] = popen(cmd, "w"); if ((uvfp[0] == NULL) | (uvfp[1] == NULL)) { fprintf(stderr, "%s: cannot open pipes to uv output\n", @@ -187,53 +191,125 @@ eval_isotropic(char *funame) } if (funame != NULL) /* need to assign Dx, Dy, Dz? */ assignD = (fundefined(funame) < 6); +#if (NSSAMP > 0) + rtrip = (char *)malloc(sqres); /* track sample triggerings */ + ryval = (float *)calloc(sqres, sizeof(float)); +#endif /* run through directions */ for (ix = 0; ix < sqres/2; ix++) { - RBFNODE *rbf = NULL; - iovec[0] = 2.*(ix+.5)/sqres - 1.; - iovec[1] = .0; - iovec[2] = input_orient * sqrt(1. - iovec[0]*iovec[0]); + const int zipsgn = (ix & 1)*2 - 1; + RBFNODE *rbf = NULL; + iovec[0] = 2.*sqfact*(ix+.5) - 1.; + iovec[1] = zipsgn*sqfact*.5; + iovec[2] = input_orient * sqrt(1. - iovec[0]*iovec[0] + - iovec[1]*iovec[1]); if (funame == NULL) rbf = advect_rbf(iovec, lobe_lim); for (ox = 0; ox < sqres; ox++) { - float last_bsdf = -1; + SDValue sdv_next; + SDsquare2disk(iovec+3, (ox+.5)*sqfact, .5*sqfact); + iovec[5] = output_orient * + sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]); + if (funame == NULL) { + eval_rbfcol(&sdv_next, rbf, iovec+3); + } else { + sdv_next.spec = c_dfcolor; + if (assignD) { + varset("Dx", '=', -iovec[3]); + varset("Dy", '=', -iovec[4]); + varset("Dz", '=', -iovec[5]); + ++eclock; + } + sdv_next.cieY = funvalue(funame, 6, iovec); + } + /* + * Super-sample when we detect a difference from before + * or after in this row, or the previous row in this column. + * Also super-sample if the previous neighbors performed + * super-sampling, regardless of any differences. + */ + memset(rtrip, 0, sqres); for (oy = 0; oy < sqres; oy++) { - SDsquare2disk(iovec+3, (ox+.5)/sqres, (oy+.5)/sqres); - iovec[5] = output_orient * + int trip; + bsdf = sdv_next.cieY; /* keeping one step ahead... */ + if (oy < sqres-1) { + SDsquare2disk(iovec+3, (ox+.5)*sqfact, (oy+1.5)*sqfact); + iovec[5] = output_orient * sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]); + } if (funame == NULL) { - SDValue sdv; - eval_rbfcol(&sdv, rbf, iovec+3); - bsdf = sdv.cieY; + SDValue sdv = sdv_next; + if (oy < sqres-1) + eval_rbfcol(&sdv_next, rbf, iovec+3); +#if (NSSAMP > 0) + trip = (abs_diff(bsdf, sdv_next.cieY) > ssamp_thresh || + (ox && abs_diff(bsdf, ryval[oy]) > ssamp_thresh) || + (oy && abs_diff(bsdf, ryval[oy-1]) > ssamp_thresh)); + if (trip | rtrip[oy] || (oy && rtrip[oy-1])) { + int ssi; + double ssa[2], sum = 0, usum = 0, vsum = 0; + /* super-sample voxel */ + for (ssi = NSSAMP; ssi--; ) { + SDmultiSamp(ssa, 2, (ssi+frandom()) * + (1./NSSAMP)); + SDsquare2disk(iovec+3, (ox+ssa[0])*sqfact, + (oy+ssa[1])*sqfact); + iovec[5] = output_orient * + sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]); + eval_rbfcol(&sdv, rbf, iovec+3); + sum += sdv.cieY; + if (rbf_colorimetry == RBCtristimulus) { + sdv.cieY /= + -2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.; + usum += 4.*sdv.spec.cx * sdv.cieY; + vsum += 9.*sdv.spec.cy * sdv.cieY; + } + } + bsdf = sum * (1./NSSAMP); + if (rbf_colorimetry == RBCtristimulus) { + uv[0] = usum / (sum+FTINY); + uv[1] = vsum / (sum+FTINY); + } + } else +#endif if (rbf_colorimetry == RBCtristimulus) { - c_ccvt(&sdv.spec, C_CSXY); uv[0] = uv[1] = 1. / (-2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.); uv[0] *= 4.*sdv.spec.cx; uv[1] *= 9.*sdv.spec.cy; } } else { - if (assignD) { - varset("Dx", '=', -iovec[3]); - varset("Dy", '=', -iovec[4]); - varset("Dz", '=', -iovec[5]); - ++eclock; + if (oy < sqres-1) { + if (assignD) { + varset("Dx", '=', -iovec[3]); + varset("Dy", '=', -iovec[4]); + varset("Dz", '=', -iovec[5]); + ++eclock; + } + sdv_next.cieY = funvalue(funame, 6, iovec); } - bsdf = funvalue(funame, 6, iovec); #if (NSSAMP > 0) - if (abs_diff(bsdf, last_bsdf) > ssamp_thresh) { + trip = (abs_diff(bsdf, sdv_next.cieY) > ssamp_thresh || + (ox && abs_diff(bsdf, ryval[oy]) > ssamp_thresh) || + (oy && abs_diff(bsdf, ryval[oy-1]) > ssamp_thresh)); + if (trip | rtrip[oy] || (oy && rtrip[oy-1])) { int ssi; - double ssa[3], ssvec[6], sum = 0; + double ssa[4], ssvec[6], sum = 0; /* super-sample voxel */ for (ssi = NSSAMP; ssi--; ) { - SDmultiSamp(ssa, 3, (ssi+frandom()) * + SDmultiSamp(ssa, 4, (ssi+frandom()) * (1./NSSAMP)); - ssvec[0] = 2.*(ix+ssa[0])/sqres - 1.; - ssvec[1] = .0; - ssvec[2] = input_orient * - sqrt(1. - ssvec[0]*ssvec[0]); - SDsquare2disk(ssvec+3, (ox+ssa[1])/sqres, - (oy+ssa[2])/sqres); + ssvec[0] = 2.*sqfact*(ix+ssa[0]) - 1.; + ssvec[1] = zipsgn*sqfact*ssa[1]; + ssvec[2] = 1. - ssvec[0]*ssvec[0] + - ssvec[1]*ssvec[1]; + if (ssvec[2] < .0) { + ssvec[1] = 0; + ssvec[2] = 1. - ssvec[0]*ssvec[0]; + } + ssvec[2] = input_orient * sqrt(ssvec[2]); + SDsquare2disk(ssvec+3, (ox+ssa[2])*sqfact, + (oy+ssa[3])*sqfact); ssvec[5] = output_orient * sqrt(1. - ssvec[3]*ssvec[3] - ssvec[4]*ssvec[4]); @@ -245,32 +321,39 @@ eval_isotropic(char *funame) } sum += funvalue(funame, 6, ssvec); } - bsdf = sum/NSSAMP; + bsdf = sum * (1./NSSAMP); } #endif } if (pctcull >= 0) - fwrite(&bsdf, sizeof(bsdf), 1, ofp); + putbinary(&bsdf, sizeof(bsdf), 1, ofp); else fprintf(ofp, "\t%.3e\n", bsdf); if (rbf_colorimetry == RBCtristimulus) { if (pctcull >= 0) { - fwrite(&uv[0], sizeof(*uv), 1, uvfp[0]); - fwrite(&uv[1], sizeof(*uv), 1, uvfp[1]); + putbinary(&uv[0], sizeof(*uv), 1, uvfp[0]); + putbinary(&uv[1], sizeof(*uv), 1, uvfp[1]); } else { fprintf(uvfp[0], "\t%.3e\n", uv[0]); fprintf(uvfp[1], "\t%.3e\n", uv[1]); } } - last_bsdf = bsdf; + if (ryval != NULL) { + rtrip[oy] = trip; + ryval[oy] = bsdf; + } } } if (rbf != NULL) free(rbf); - prog_show((ix+1.)*(2./sqres)); + prog_show((ix+1.)*(2.*sqfact)); } prog_done(); + if (ryval != NULL) { + free(rtrip); + free(ryval); + } if (pctcull >= 0) { /* finish output */ if (pclose(ofp)) { fprintf(stderr, "%s: error running rttree_reduce on Y\n", @@ -313,6 +396,9 @@ static void eval_anisotropic(char *funame) { const int sqres = 1<= 0) { const char *avgopt = (input_orient>0 ^ output_orient>0) - ? "" : " -a"; + ? "" : recip; sprintf(cmd, "rttree_reduce%s -h -ff -r 4 -t %f -g %d > %s", avgopt, pctcull, samp_order, create_component_file(0)); @@ -378,55 +464,120 @@ eval_anisotropic(char *funame) } if (funame != NULL) /* need to assign Dx, Dy, Dz? */ assignD = (fundefined(funame) < 6); +#if (NSSAMP > 0) + rtrip = (char *)malloc(sqres); /* track sample triggerings */ + ryval = (float *)calloc(sqres, sizeof(float)); +#endif /* run through directions */ for (ix = 0; ix < sqres; ix++) for (iy = 0; iy < sqres; iy++) { RBFNODE *rbf = NULL; /* Klems reversal */ - SDsquare2disk(iovec, 1.-(ix+.5)/sqres, 1.-(iy+.5)/sqres); + SDsquare2disk(iovec, 1.-(ix+.5)*sqfact, 1.-(iy+.5)*sqfact); iovec[2] = input_orient * sqrt(1. - iovec[0]*iovec[0] - iovec[1]*iovec[1]); if (funame == NULL) rbf = advect_rbf(iovec, lobe_lim); for (ox = 0; ox < sqres; ox++) { - float last_bsdf = -1; + SDValue sdv_next; + SDsquare2disk(iovec+3, (ox+.5)*sqfact, .5*sqfact); + iovec[5] = output_orient * + sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]); + if (funame == NULL) { + eval_rbfcol(&sdv_next, rbf, iovec+3); + } else { + sdv_next.spec = c_dfcolor; + if (assignD) { + varset("Dx", '=', -iovec[3]); + varset("Dy", '=', -iovec[4]); + varset("Dz", '=', -iovec[5]); + ++eclock; + } + sdv_next.cieY = funvalue(funame, 6, iovec); + } + /* + * Super-sample when we detect a difference from before + * or after in this row, or the previous row in this column. + * Also super-sample if the previous neighbors performed + * super-sampling, regardless of any differences. + */ + memset(rtrip, 0, sqres); for (oy = 0; oy < sqres; oy++) { - SDsquare2disk(iovec+3, (ox+.5)/sqres, (oy+.5)/sqres); - iovec[5] = output_orient * + int trip; + bsdf = sdv_next.cieY; /* keeping one step ahead... */ + if (oy < sqres-1) { + SDsquare2disk(iovec+3, (ox+.5)*sqfact, (oy+1.5)*sqfact); + iovec[5] = output_orient * sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]); + } if (funame == NULL) { - SDValue sdv; - eval_rbfcol(&sdv, rbf, iovec+3); - bsdf = sdv.cieY; + SDValue sdv = sdv_next; + if (oy < sqres-1) + eval_rbfcol(&sdv_next, rbf, iovec+3); +#if (NSSAMP > 0) + trip = (abs_diff(bsdf, sdv_next.cieY) > ssamp_thresh || + (ox && abs_diff(bsdf, ryval[oy]) > ssamp_thresh) || + (oy && abs_diff(bsdf, ryval[oy-1]) > ssamp_thresh)); + if (trip | rtrip[oy] || (oy && rtrip[oy-1])) { + int ssi; + double ssa[2], sum = 0, usum = 0, vsum = 0; + /* super-sample voxel */ + for (ssi = NSSAMP; ssi--; ) { + SDmultiSamp(ssa, 2, (ssi+frandom()) * + (1./NSSAMP)); + SDsquare2disk(iovec+3, (ox+ssa[0])*sqfact, + (oy+ssa[1])*sqfact); + iovec[5] = output_orient * + sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]); + eval_rbfcol(&sdv, rbf, iovec+3); + sum += sdv.cieY; + if (rbf_colorimetry == RBCtristimulus) { + sdv.cieY /= + -2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.; + usum += 4.*sdv.spec.cx * sdv.cieY; + vsum += 9.*sdv.spec.cy * sdv.cieY; + } + } + bsdf = sum * (1./NSSAMP); + if (rbf_colorimetry == RBCtristimulus) { + uv[0] = usum / (sum+FTINY); + uv[1] = vsum / (sum+FTINY); + } + } else +#endif if (rbf_colorimetry == RBCtristimulus) { - c_ccvt(&sdv.spec, C_CSXY); uv[0] = uv[1] = 1. / (-2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.); uv[0] *= 4.*sdv.spec.cx; uv[1] *= 9.*sdv.spec.cy; } } else { - if (assignD) { - varset("Dx", '=', -iovec[3]); - varset("Dy", '=', -iovec[4]); - varset("Dz", '=', -iovec[5]); - ++eclock; + if (oy < sqres-1) { + if (assignD) { + varset("Dx", '=', -iovec[3]); + varset("Dy", '=', -iovec[4]); + varset("Dz", '=', -iovec[5]); + ++eclock; + } + sdv_next.cieY = funvalue(funame, 6, iovec); } - bsdf = funvalue(funame, 6, iovec); #if (NSSAMP > 0) - if (abs_diff(bsdf, last_bsdf) > ssamp_thresh) { + trip = (abs_diff(bsdf, sdv_next.cieY) > ssamp_thresh || + (ox && abs_diff(bsdf, ryval[oy]) > ssamp_thresh) || + (oy && abs_diff(bsdf, ryval[oy-1]) > ssamp_thresh)); + if (trip | rtrip[oy] || (oy && rtrip[oy-1])) { int ssi; double ssa[4], ssvec[6], sum = 0; /* super-sample voxel */ for (ssi = NSSAMP; ssi--; ) { SDmultiSamp(ssa, 4, (ssi+frandom()) * (1./NSSAMP)); - SDsquare2disk(ssvec, 1.-(ix+ssa[0])/sqres, - 1.-(iy+ssa[1])/sqres); + SDsquare2disk(ssvec, 1.-(ix+ssa[0])*sqfact, + 1.-(iy+ssa[1])*sqfact); ssvec[2] = input_orient * sqrt(1. - ssvec[0]*ssvec[0] - ssvec[1]*ssvec[1]); - SDsquare2disk(ssvec+3, (ox+ssa[2])/sqres, - (oy+ssa[3])/sqres); + SDsquare2disk(ssvec+3, (ox+ssa[2])*sqfact, + (oy+ssa[3])*sqfact); ssvec[5] = output_orient * sqrt(1. - ssvec[3]*ssvec[3] - ssvec[4]*ssvec[4]); @@ -438,32 +589,39 @@ eval_anisotropic(char *funame) } sum += funvalue(funame, 6, ssvec); } - bsdf = sum/NSSAMP; + bsdf = sum * (1./NSSAMP); } #endif } if (pctcull >= 0) - fwrite(&bsdf, sizeof(bsdf), 1, ofp); + putbinary(&bsdf, sizeof(bsdf), 1, ofp); else fprintf(ofp, "\t%.3e\n", bsdf); if (rbf_colorimetry == RBCtristimulus) { if (pctcull >= 0) { - fwrite(&uv[0], sizeof(*uv), 1, uvfp[0]); - fwrite(&uv[1], sizeof(*uv), 1, uvfp[1]); + putbinary(&uv[0], sizeof(*uv), 1, uvfp[0]); + putbinary(&uv[1], sizeof(*uv), 1, uvfp[1]); } else { fprintf(uvfp[0], "\t%.3e\n", uv[0]); fprintf(uvfp[1], "\t%.3e\n", uv[1]); } } - last_bsdf = bsdf; + if (ryval != NULL) { + rtrip[oy] = trip; + ryval[oy] = bsdf; + } } - } if (rbf != NULL) free(rbf); prog_show((ix*sqres+iy+1.)/(sqres*sqres)); } + } prog_done(); + if (ryval != NULL) { + free(rtrip); + free(ryval); + } if (pctcull >= 0) { /* finish output */ if (pclose(ofp)) { fprintf(stderr, "%s: error running rttree_reduce on Y\n", @@ -565,21 +723,37 @@ main(int argc, char *argv[]) esupport &= ~(E_INCHAN|E_OUTCHAN); scompile("PI:3.14159265358979323846", NULL, 0); biggerlib(); - for (i = 1; i < argc-1 && (argv[i][0] == '-') | (argv[i][0] == '+'); i++) + for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++) switch (argv[i][1]) { /* get options */ case 'e': scompile(argv[++i], NULL, 0); + if (single_plane_incident < 0) + single_plane_incident = 0; break; case 'f': if (!argv[i][2]) { if (strchr(argv[++i], '=') != NULL) { add_wbsdf("-f", 1); add_wbsdf(argv[i], 1); - } else - fcompile(argv[i]); + } else { + 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); + if (single_plane_incident < 0) + single_plane_incident = 0; + } } else dofwd = (argv[i][0] == '+'); break; + case 'a': + recip = (argv[i][0] == '+') ? " -a" : ""; + break; case 'b': dobwd = (argv[i][0] == '+'); break; @@ -633,6 +807,7 @@ main(int argc, char *argv[]) goto userr; } ++eclock; + ssamp_thresh *= 0.5; /* lower sampling threshold */ add_wbsdf("-a", 1); add_wbsdf(tfmt[single_plane_incident], 1); if (dofwd) { @@ -695,10 +870,10 @@ main(int argc, char *argv[]) return(wrap_up()); userr: fprintf(stderr, - "Usage: %s [-g Nlog2][-t pctcull][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", + "Usage: %s [{+|-}a][-g Nlog2][-t pctcull][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname); fprintf(stderr, - " or: %s -t{3|4} [-g Nlog2][-t pctcull][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n", + " or: %s -t{3|4} [{+|-}a][-g Nlog2][-t pctcull][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n", progname); return(1); }