--- ray/src/cv/pabopto2bsdf.c 2014/03/20 16:54:23 2.20 +++ ray/src/cv/pabopto2bsdf.c 2021/04/05 19:56:18 2.40 @@ -1,35 +1,39 @@ #ifndef lint -static const char RCSid[] = "$Id: pabopto2bsdf.c,v 2.20 2014/03/20 16:54:23 greg Exp $"; +static const char RCSid[] = "$Id: pabopto2bsdf.c,v 2.40 2021/04/05 19:56:18 greg Exp $"; #endif /* * Load measured BSDF data in PAB-Opto format. + * Assumes that surface-normal (Z-axis) faces into room unless -t option given. * * G. Ward */ #define _USE_MATH_DEFINES -#include #include -#include #include #include +#include "rtio.h" #include "platform.h" #include "bsdfrep.h" -#include "resolu.h" /* global argv[0] */ char *progname; typedef struct { const char *fname; /* input file path */ - double theta, phi; /* input angles */ + double theta, phi; /* input angles (degrees) */ + double up_phi; /* azimuth for "up" direction */ int igp[2]; /* input grid position */ int isDSF; /* data is DSF (rather than BSDF)? */ + int nspec; /* number of spectral samples */ long dstart; /* data start offset in file */ } PGINPUT; PGINPUT *inpfile; /* input files sorted by incidence */ -int ninpfiles; /* number of input files */ +int rev_orient = 0; /* shall we reverse surface orientation? */ + +double lim_graze = 0; /* limit scattering near grazing (deg above) */ + /* Compare incident angles */ static int cmp_indir(const void *p1, const void *p2) @@ -44,12 +48,28 @@ cmp_indir(const void *p1, const void *p2) return(inp1->igp[0] - inp2->igp[0]); } +/* Assign grid position from theta and phi */ +static void +set_grid_pos(PGINPUT *pip) +{ + FVECT dv; + + if (pip->theta <= FTINY) { + pip->igp[0] = pip->igp[1] = grid_res/2 - 1; + return; + } + dv[2] = sin(M_PI/180.*pip->theta); + dv[0] = cos(M_PI/180.*pip->phi)*dv[2]; + dv[1] = sin(M_PI/180.*pip->phi)*dv[2]; + dv[2] = sqrt(1. - dv[2]*dv[2]); + pos_from_vec(pip->igp, dv); +} + /* Prepare a PAB-Opto input file by reading its header */ static int init_pabopto_inp(const int i, const char *fname) { FILE *fp = fopen(fname, "r"); - FVECT dv; char buf[2048]; int c; @@ -60,10 +80,12 @@ init_pabopto_inp(const int i, const char *fname) } inpfile[i].fname = fname; inpfile[i].isDSF = -1; + inpfile[i].nspec = 0; + inpfile[i].up_phi = 0; inpfile[i].theta = inpfile[i].phi = -10001.; /* read header information */ while ((c = getc(fp)) == '#' || c == EOF) { - char typ[32]; + char typ[64]; if (fgets(buf, sizeof(buf), fp) == NULL) { fputs(fname, stderr); fputs(": unexpected EOF\n", stderr); @@ -72,16 +94,24 @@ init_pabopto_inp(const int i, const char *fname) } if (sscanf(buf, "sample_name \"%[^\"]\"", bsdf_name) == 1) continue; + if (sscanf(buf, "colorimetry: %s", typ) == 1) { + if (!strcasecmp(typ, "CIE-XYZ")) + inpfile[i].nspec = 3; + else if (!strcasecmp(typ, "CIE-Y")) + inpfile[i].nspec = 1; + continue; + } if (sscanf(buf, "format: theta phi %s", typ) == 1) { - if (!strcasecmp(typ, "DSF")) { + if (!strcasecmp(typ, "DSF")) inpfile[i].isDSF = 1; - continue; - } - if (!strcasecmp(typ, "BSDF")) { + else if (!strcasecmp(typ, "BSDF") || + !strcasecmp(typ, "BRDF") || + !strcasecmp(typ, "BTDF")) inpfile[i].isDSF = 0; - continue; - } + continue; } + if (sscanf(buf, "upphi %lf", &inpfile[i].up_phi) == 1) + continue; if (sscanf(buf, "intheta %lf", &inpfile[i].theta) == 1) continue; if (sscanf(buf, "inphi %lf", &inpfile[i].phi) == 1) @@ -102,12 +132,14 @@ init_pabopto_inp(const int i, const char *fname) fputs(": unknown incident angle\n", stderr); return(0); } + if (rev_orient) { /* reverse Z-axis to face outside */ + inpfile[i].theta = 180. - inpfile[i].theta; + inpfile[i].phi = 360. - inpfile[i].phi; + } + /* convert to Y-up orientation */ + inpfile[i].phi += 90.-inpfile[i].up_phi; /* convert angle to grid position */ - dv[2] = sin(M_PI/180.*inpfile[i].theta); - dv[0] = cos(M_PI/180.*inpfile[i].phi)*dv[2]; - dv[1] = sin(M_PI/180.*inpfile[i].phi)*dv[2]; - dv[2] = sqrt(1. - dv[2]*dv[2]); - pos_from_vec(inpfile[i].igp, dv); + set_grid_pos(&inpfile[i]); return(1); } @@ -116,7 +148,7 @@ static int add_pabopto_inp(const int i) { FILE *fp = fopen(inpfile[i].fname, "r"); - double theta_out, phi_out, val; + double theta_out, phi_out, val[3]; int n, c; if (fp == NULL || fseek(fp, inpfile[i].dstart, 0) == EOF) { @@ -129,17 +161,35 @@ add_pabopto_inp(const int i) if (i) /* process previous incidence */ make_rbfrep(); #ifdef DEBUG - fprintf(stderr, "New incident (theta,phi)=(%f,%f)\n", + fprintf(stderr, "New incident (theta,phi)=(%.1f,%.1f)\n", inpfile[i].theta, inpfile[i].phi); #endif + if (inpfile[i].nspec) + set_spectral_samples(inpfile[i].nspec); new_bsdf_data(inpfile[i].theta, inpfile[i].phi); } #ifdef DEBUG fprintf(stderr, "Loading measurements from '%s'...\n", inpfile[i].fname); #endif /* read scattering data */ - while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) + while (fscanf(fp, "%lf %lf %lf", &theta_out, &phi_out, val) == 3) { + for (n = 1; n < inpfile[i].nspec; n++) + if (fscanf(fp, "%lf", val+n) != 1) { + fprintf(stderr, "%s: warning: unexpected EOF\n", + inpfile[i].fname); + fclose(fp); + return(1); + } + /* check if scatter angle is too low */ + if (fabs(theta_out - 90.) < lim_graze-FTINY) + continue; + if (rev_orient) { /* reverse Z-axis to face outside? */ + theta_out = 180. - theta_out; + phi_out = 360. - phi_out; + } + phi_out += 90.-inpfile[i].up_phi; add_bsdf_data(theta_out, phi_out, val, inpfile[i].isDSF); + } n = 0; while ((c = getc(fp)) != EOF) n += !isspace(c); @@ -152,53 +202,136 @@ add_pabopto_inp(const int i) } #ifndef TEST_MAIN + +#define SYM_ILL '?' /* illegal symmetry value */ +#define SYM_ISO 'I' /* isotropic */ +#define SYM_QUAD 'Q' /* quadrilateral symmetry */ +#define SYM_BILAT 'B' /* bilateral symmetry */ +#define SYM_ANISO 'A' /* anisotropic */ +#define SYM_UP 'U' /* "up-down" (180°) symmetry */ + +static const char quadrant_rep[16][16] = { + "in-plane","0-90","90-180","0-180", + "180-270","0-90+180-270","90-270", + "0-270","270-360","270-90", + "90-180+270-360","270-180","180-360", + "180-90","90-360","0-360" + }; +static const char quadrant_sym[16] = { + SYM_ISO, SYM_QUAD, SYM_QUAD, SYM_BILAT, + SYM_QUAD, SYM_ILL, SYM_BILAT, SYM_ILL, + SYM_QUAD, SYM_BILAT, SYM_ILL, SYM_ILL, + SYM_BILAT, SYM_ILL, SYM_ILL, SYM_ANISO + }; + /* Read in PAB-Opto BSDF files and output RBF interpolant */ int main(int argc, char *argv[]) { extern int nprocs; - int i; - /* start header */ - SET_FILE_BINARY(stdout); - newheader("RADIANCE", stdout); - printargs(argc, argv, stdout); - fputnow(stdout); + int auto_grazing = 0; + const char *symmetry = "0"; + int ninpfiles, totinc; + int a, i; + progname = argv[0]; /* get options */ - while (argc > 2 && argv[1][0] == '-') { - switch (argv[1][1]) { + for (a = 1; a < argc && argv[a][0] == '-'; a++) + switch (argv[a][1]) { + case 't': + rev_orient = !rev_orient; + break; case 'n': - nprocs = atoi(argv[2]); + nprocs = atoi(argv[++a]); break; + case 's': + symmetry = argv[++a]; + break; + case 'g': + if (toupper(argv[a+1][0]) == 'A') + auto_grazing = 1; + else + lim_graze = atof(argv[a+1]); + ++a; + break; default: goto userr; } - argv += 2; argc -= 2; - } - /* initialize & sort inputs */ - ninpfiles = argc - 1; + totinc = ninpfiles = argc - a; /* initialize & sort inputs */ if (ninpfiles < 2) goto userr; - inpfile = (PGINPUT *)malloc(sizeof(PGINPUT)*ninpfiles); + if (toupper(symmetry[0]) == SYM_UP) /* special case for "up" symmetry */ + totinc += ninpfiles; + inpfile = (PGINPUT *)malloc(sizeof(PGINPUT)*totinc); if (inpfile == NULL) return(1); - for (i = 0; i < ninpfiles; i++) - if (!init_pabopto_inp(i, argv[i+1])) + if (auto_grazing) + lim_graze = 90.; + for (i = 0; i < ninpfiles; i++) { + if (!init_pabopto_inp(i, argv[a+i])) return(1); - qsort(inpfile, ninpfiles, sizeof(PGINPUT), &cmp_indir); + if (auto_grazing && fabs(inpfile[i].theta - 90.) < lim_graze) + lim_graze = fabs(inpfile[i].theta - 90.); + } + for (i = ninpfiles; i < totinc; i++) { /* copy for "up" symmetry */ + inpfile[i] = inpfile[i-ninpfiles]; + inpfile[i].phi += 180.; /* invert duplicate data */ + inpfile[i].up_phi -= 180.; + set_grid_pos(&inpfile[i]); /* grid location for sorting */ + } + qsort(inpfile, totinc, sizeof(PGINPUT), cmp_indir); /* compile measurements */ - for (i = 0; i < ninpfiles; i++) + for (i = 0; i < totinc; i++) if (!add_pabopto_inp(i)) return(1); make_rbfrep(); /* process last data set */ + /* check input symmetry */ + switch (toupper(symmetry[0])) { + case '0': /* unspecified symmetry */ + if (quadrant_sym[inp_coverage] != SYM_ILL) + break; /* anything legal goes */ + fprintf(stderr, "%s: unsupported phi coverage (%s)\n", + progname, quadrant_rep[inp_coverage]); + return(1); + case SYM_UP: /* faux "up" symmetry */ + if (quadrant_sym[inp_coverage] == SYM_ANISO) + break; + /* fall through */ + case SYM_ISO: /* usual symmetry types */ + case SYM_QUAD: + case SYM_BILAT: + case SYM_ANISO: + if (quadrant_sym[inp_coverage] == toupper(symmetry[0])) + break; /* matches spec */ + fprintf(stderr, + "%s: phi coverage (%s) does not match requested '%s' symmetry\n", + progname, quadrant_rep[inp_coverage], symmetry); + return(1); + default: + fprintf(stderr, + "%s: -s option must be Isotropic, Quadrilateral, Bilateral, Up, or Anisotropic\n", + progname); + return(1); + } +#ifdef DEBUG + fprintf(stderr, "Input phi coverage (%s) has '%c' symmetry\n", + quadrant_rep[inp_coverage], + quadrant_sym[inp_coverage]); +#endif build_mesh(); /* create interpolation */ - save_bsdf_rep(stdout); /* write it out */ + SET_FILE_BINARY(stdout); /* start header */ + newheader("RADIANCE", stdout); + printargs(argc, argv, stdout); + fputnow(stdout); + save_bsdf_rep(stdout); /* complete header + data */ return(0); userr: - fprintf(stderr, "Usage: %s [-n nproc] meas1.dat meas2.dat .. > bsdf.sir\n", + fprintf(stderr, "Usage: %s [-t][-n nproc][-s symmetry][-g angle|'A'] meas1.dat meas2.dat .. > bsdf.sir\n", progname); return(1); } -#else + +#else /* TEST_MAIN */ + /* Test main produces a Radiance model from the given input file */ int main(int argc, char *argv[]) @@ -210,16 +343,19 @@ main(int argc, char *argv[]) FVECT dir; int i, j, n; + progname = argv[0]; if (argc != 2) { - fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]); + fprintf(stderr, "Usage: %s input.dat > output.rad\n", progname); return(1); } - ninpfiles = 1; inpfile = &pginp; if (!init_pabopto_inp(0, argv[1]) || !add_pabopto_inp(0)) return(1); /* reduce data set */ - make_rbfrep(); + if (make_rbfrep() == NULL) { + fprintf(stderr, "%s: nothing to plot!\n", progname); + exit(1); + } #ifdef DEBUG fprintf(stderr, "Minimum BSDF = %.4f\n", bsdf_min); #endif @@ -227,12 +363,12 @@ main(int argc, char *argv[]) #if 1 /* produce spheres at meas. */ puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n"); n = 0; - for (i = 0; i < GRIDRES; i++) - for (j = 0; j < GRIDRES; j++) + for (i = 0; i < grid_res; i++) + for (j = 0; j < grid_res; j++) if (dsf_grid[i][j].sum.n > 0) { ovec_from_pos(dir, i, j); bsdf = dsf_grid[i][j].sum.v / - (dsf_grid[i][j].sum.n*output_orient*dir[2]); + ((double)dsf_grid[i][j].sum.n*output_orient*dir[2]); if (bsdf <= bsdf_min*.6) continue; bsdf = log(bsdf + 1e-5) - min_log; @@ -248,7 +384,7 @@ main(int argc, char *argv[]) for (n = 0; n < dsf_list->nrbf; n++) { RBFVAL *rbf = &dsf_list->rbfa[n]; ovec_from_pos(dir, rbf->gx, rbf->gy); - bsdf = eval_rbfrep(dsf_list, dir) / (output_orient*dir[2]); + bsdf = eval_rbfrep(dsf_list, dir); bsdf = log(bsdf + 1e-5) - min_log; printf("red sphere p%04d\n0\n0\n", ++n); printf("4 %.6g %.6g %.6g %.6g\n\n", @@ -259,16 +395,16 @@ main(int argc, char *argv[]) #if 1 /* output continuous surface */ puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 1\n"); fflush(stdout); - sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1); + sprintf(buf, "gensurf tgreen bsdf - - - %d %d", grid_res-1, grid_res-1); pfp = popen(buf, "w"); if (pfp == NULL) { - fprintf(stderr, "%s: cannot open '| %s'\n", argv[0], buf); + fprintf(stderr, "%s: cannot open '| %s'\n", progname, buf); return(1); } - for (i = 0; i < GRIDRES; i++) - for (j = 0; j < GRIDRES; j++) { + for (i = 0; i < grid_res; i++) + for (j = 0; j < grid_res; j++) { ovec_from_pos(dir, i, j); - bsdf = eval_rbfrep(dsf_list, dir) / (output_orient*dir[2]); + bsdf = eval_rbfrep(dsf_list, dir); bsdf = log(bsdf + 1e-5) - min_log; fprintf(pfp, "%.8e %.8e %.8e\n", dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf); @@ -278,4 +414,5 @@ main(int argc, char *argv[]) #endif return(0); } -#endif + +#endif /* TEST_MAIN */