--- ray/src/cv/bsdf2klems.c 2013/04/21 23:01:14 2.1 +++ ray/src/cv/bsdf2klems.c 2013/08/02 20:56:19 2.8 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf2klems.c,v 2.1 2013/04/21 23:01:14 greg Exp $"; +static const char RCSid[] = "$Id: bsdf2klems.c,v 2.8 2013/08/02 20:56:19 greg Exp $"; #endif /* * Load measured BSDF interpolant and write out as XML file with Klems matrix. @@ -17,6 +17,8 @@ static const char RCSid[] = "$Id: bsdf2klems.c,v 2.1 2 #include "calcomp.h" #include "bsdfrep.h" #include "bsdf_m.h" + /* assumed maximum # Klems patches */ +#define MAXPATCHES 145 /* global argv[0] */ char *progname; /* selected basis function name */ @@ -36,17 +38,10 @@ get_basis(const char *bn) return NULL; } -/* Output XML prologue to stdout */ +/* Output XML header to stdout */ static void -xml_prologue(int ac, char *av[]) +xml_header(int ac, char *av[]) { - ANGLE_BASIS *abp = get_basis(kbasis); - int i; - - if (abp == NULL) { - fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis); - exit(1); - } puts(""); puts(""); fputs(""); +} + +/* Output XML prologue to stdout */ +static void +xml_prologue(const SDData *sd) +{ + const char *matn = (sd && sd->matn[0]) ? sd->matn : "Name"; + const char *makr = (sd && sd->makr[0]) ? sd->makr : "Manufacturer"; + ANGLE_BASIS *abp = get_basis(kbasis); + int i; + + if (abp == NULL) { + fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis); + exit(1); + } puts("System"); puts("BSDF"); puts(""); puts(""); puts("\t"); - puts("\t\tName"); - puts("\t\tManufacturer"); + printf("\t\t%s\n", matn); + printf("\t\t%s\n", makr); + if (sd && sd->dim[2] > .001) { + printf("\t\t%.3f\n", sd->dim[2]); + printf("\t\t%.3f\n", sd->dim[0]); + printf("\t\t%.3f\n", sd->dim[1]); + } puts("\t\tOther"); puts("\t"); + if (sd && sd->mgf != NULL) { + puts("\t"); + puts("\t\t"); + fputs(sd->mgf, stdout); + puts(""); + puts("\t"); + } puts("\t"); puts("\t\tColumns"); puts("\t\t"); @@ -73,7 +95,7 @@ xml_prologue(int ac, char *av[]) printf("\t\t\t%g\n", i ? .5*(abp->lat[i].tmin + abp->lat[i+1].tmin) : .0 ); - printf("\t\t\t%d", abp->lat[i].nphis); + printf("\t\t\t%d\n", abp->lat[i].nphis); puts("\t\t\t"); printf("\t\t\t\t%g\n", abp->lat[i].tmin); printf("\t\t\t\t%g\n", abp->lat[i+1].tmin); @@ -127,11 +149,12 @@ xml_epilogue(void) puts(""); } -/* Load and resample XML BSDF description */ +/* Load and resample XML BSDF description using Klems basis */ static void eval_bsdf(const char *fname) { ANGLE_BASIS *abp = get_basis(kbasis); + float *trans_mtx = NULL; SDData bsd; SDError ec; FVECT vin, vout; @@ -142,6 +165,7 @@ eval_bsdf(const char *fname) SDclearBSDF(&bsd, fname); /* load BSDF file */ if ((ec = SDloadFile(&bsd, fname)) != SDEnone) goto err; + xml_prologue(&bsd); /* pass geometry */ /* front reflection */ if (bsd.rf != NULL || bsd.rLambFront.cieY > .002) { input_orient = 1; output_orient = 1; @@ -151,7 +175,7 @@ eval_bsdf(const char *fname) sum = 0; /* average over patches */ for (n = npsamps; n-- > 0; ) { fo_getvec(vout, j+(n+frandom())/npsamps, abp); - fi_getvec(vin, i+(n+frandom())/npsamps, abp); + fi_getvec(vin, i+urand(n), abp); ec = SDevalBSDF(&sv, vout, vin, &bsd); if (ec != SDEnone) goto err; @@ -172,7 +196,7 @@ eval_bsdf(const char *fname) sum = 0; /* average over patches */ for (n = npsamps; n-- > 0; ) { bo_getvec(vout, j+(n+frandom())/npsamps, abp); - bi_getvec(vin, i+(n+frandom())/npsamps, abp); + bi_getvec(vin, i+urand(n), abp); ec = SDevalBSDF(&sv, vout, vin, &bsd); if (ec != SDEnone) goto err; @@ -186,6 +210,9 @@ eval_bsdf(const char *fname) } /* front transmission */ if (bsd.tf != NULL || bsd.tLamb.cieY > .002) { + if (bsd.tb == NULL) + trans_mtx = (float *)malloc(sizeof(float) * + abp->nangles*abp->nangles); input_orient = 1; output_orient = -1; data_prologue(); for (j = 0; j < abp->nangles; j++) { @@ -193,55 +220,80 @@ eval_bsdf(const char *fname) sum = 0; /* average over patches */ for (n = npsamps; n-- > 0; ) { bo_getvec(vout, j+(n+frandom())/npsamps, abp); - fi_getvec(vin, i+(n+frandom())/npsamps, abp); + fi_getvec(vin, i+urand(n), abp); ec = SDevalBSDF(&sv, vout, vin, &bsd); if (ec != SDEnone) goto err; sum += sv.cieY; } printf("\t%.3e\n", sum/npsamps); + if (trans_mtx != NULL) + trans_mtx[j*abp->nangles + i] = sum/npsamps; } putchar('\n'); /* extra space between rows */ } data_epilogue(); } /* back transmission */ - if (bsd.tb != NULL) { + if (bsd.tb != NULL || trans_mtx != NULL) { + if (bsd.tf == NULL) + trans_mtx = (float *)malloc(sizeof(float) * + abp->nangles*abp->nangles); input_orient = -1; output_orient = 1; data_prologue(); for (j = 0; j < abp->nangles; j++) { - for (i = 0; i < abp->nangles; i++) { - sum = 0; /* average over patches */ - for (n = npsamps; n-- > 0; ) { - fo_getvec(vout, j+(n+frandom())/npsamps, abp); - bi_getvec(vin, i+(n+frandom())/npsamps, abp); - ec = SDevalBSDF(&sv, vout, vin, &bsd); - if (ec != SDEnone) + for (i = 0; i < abp->nangles; i++) + if (bsd.tb != NULL) { /* use tb if we have it */ + sum = 0; /* average over patches */ + for (n = npsamps; n-- > 0; ) { + fo_getvec(vout, j+(n+frandom())/npsamps, abp); + bi_getvec(vin, i+urand(n), abp); + ec = SDevalBSDF(&sv, vout, vin, &bsd); + if (ec != SDEnone) goto err; - sum += sv.cieY; + sum += sv.cieY; + } + printf("\t%.3e\n", sum/npsamps); + if (trans_mtx != NULL) + trans_mtx[i*abp->nangles + j] = sum/npsamps; + } else { /* else transpose tf */ + printf("\t%.3e\n", trans_mtx[i*abp->nangles + j]); } - printf("\t%.3e\n", sum/npsamps); - } putchar('\n'); /* extra space between rows */ } data_epilogue(); } + /* derived front transmission */ + if (bsd.tf == NULL && trans_mtx != NULL) { + input_orient = 1; output_orient = -1; + data_prologue(); + for (j = 0; j < abp->nangles; j++) { + for (i = 0; i < abp->nangles; i++) + printf("\t%.3e\n", trans_mtx[j*abp->nangles + i]); + putchar('\n'); /* extra space between rows */ + } + data_epilogue(); + } SDfreeBSDF(&bsd); /* all done */ + if (trans_mtx != NULL) + free(trans_mtx); return; err: SDreportError(ec, stderr); exit(1); } -/* Interpolate and output a BSDF function */ +/* Interpolate and output a BSDF function using Klems basis */ static void eval_function(char *funame) { ANGLE_BASIS *abp = get_basis(kbasis); + int assignD = (fundefined(funame) < 6); double iovec[6]; double sum; int i, j, n; + initurand(npsamps); data_prologue(); /* begin output */ for (j = 0; j < abp->nangles; j++) { /* run through directions */ for (i = 0; i < abp->nangles; i++) { @@ -253,10 +305,16 @@ eval_function(char *funame) bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp); if (input_orient > 0) - fi_getvec(iovec, j+(n+frandom())/npsamps, abp); + fi_getvec(iovec, i+urand(n), abp); else - bi_getvec(iovec, j+(n+frandom())/npsamps, abp); + bi_getvec(iovec, i+urand(n), abp); + if (assignD) { + varset("Dx", '=', -iovec[3]); + varset("Dy", '=', -iovec[4]); + varset("Dz", '=', -iovec[5]); + ++eclock; + } sum += funvalue(funame, 6, iovec); } printf("\t%.3e\n", sum/npsamps); @@ -271,12 +329,47 @@ static void eval_rbf(void) { ANGLE_BASIS *abp = get_basis(kbasis); - double iovec[6]; + float bsdfarr[MAXPATCHES*MAXPATCHES]; + FVECT vin, vout; + RBFNODE *rbf; double sum; int i, j, n; + /* sanity check */ + if (abp->nangles > MAXPATCHES) { + fprintf(stderr, "%s: too many patches!\n", progname); + exit(1); + } + data_prologue(); /* begin output */ + for (i = 0; i < abp->nangles; i++) { + 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); -fprintf(stder, "%s: RBF evaluation currently unimplemented\n", progname); -exit(1); + rbf = advect_rbf(vin); /* compute radial basis func */ + + for (j = 0; j < abp->nangles; j++) { + sum = 0; /* sample over exiting patch */ + for (n = npsamps; n--; ) { + if (output_orient > 0) + fo_getvec(vout, j+(n+frandom())/npsamps, abp); + else + bo_getvec(vout, j+(n+frandom())/npsamps, abp); + + sum += eval_rbfrep(rbf, vout) / vout[2]; + } + bsdfarr[j*abp->nangles + i] = sum*output_orient/npsamps; + } + if (rbf != NULL) + free(rbf); + } + n = 0; /* write out our matrix */ + for (j = 0; j < abp->nangles; j++) { + for (i = 0; i < abp->nangles; i++) + printf("\t%.3e\n", bsdfarr[n++]); + putchar('\n'); + } + data_epilogue(); /* finish output */ } /* Read in BSDF and interpolate as Klems matrix representation */ @@ -292,7 +385,7 @@ 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 'n': npsamps = atoi(argv[++i]); @@ -327,9 +420,13 @@ main(int argc, char *argv[]) fprintf(stderr, "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n", progname); + fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n", + progname); goto userr; } - xml_prologue(argc, argv); /* start XML output */ + ++eclock; + xml_header(argc, argv); /* start XML output */ + xml_prologue(NULL); if (dofwd) { input_orient = -1; output_orient = -1; @@ -347,9 +444,10 @@ main(int argc, char *argv[]) xml_epilogue(); /* finish XML output & exit */ return(0); } - if (i == argc-1 && (cp = strstr(argv[i], ".xml")) != NULL && - strlen(cp) == 4) { /* XML input? */ - xml_prologue(argc, argv); /* start XML output */ + /* XML input? */ + if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] && + !strcasecmp(cp, ".xml")) { + xml_header(argc, argv); /* start XML output */ eval_bsdf(argv[i]); /* load & resample BSDF */ xml_epilogue(); /* finish XML output & exit */ return(0); @@ -366,8 +464,10 @@ main(int argc, char *argv[]) if (!load_bsdf_rep(fpin)) return(1); fclose(fpin); - if (!nbsdf++) /* start XML on first dist. */ - xml_prologue(argc, argv); + if (!nbsdf++) { /* start XML on first dist. */ + xml_header(argc, argv); + xml_prologue(NULL); + } eval_rbf(); } xml_epilogue(); /* finish XML output & exit */ @@ -376,17 +476,16 @@ main(int argc, char *argv[]) SET_FILE_BINARY(stdin); /* load from stdin */ if (!load_bsdf_rep(stdin)) return(1); - xml_prologue(argc, argv); /* start XML output */ + xml_header(argc, argv); /* start XML output */ + xml_prologue(NULL); eval_rbf(); /* resample dist. */ xml_epilogue(); /* finish XML output & exit */ return(0); userr: fprintf(stderr, - "Usage: %s [-n spp][-h|-q][bsdf.sir ..] > bsdf.xml\n", - progname); + "Usage: %s [-n spp][-h|-q][bsdf.sir ..] > bsdf.xml\n", progname); fprintf(stderr, - " or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", - progname); + " or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", progname); fprintf(stderr, " or: %s [-n spp][-h|-q][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n", progname);