--- ray/src/util/wrapBSDF.c 2015/02/13 20:49:59 2.3 +++ ray/src/util/wrapBSDF.c 2015/05/29 07:45:48 2.16 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: wrapBSDF.c,v 2.3 2015/02/13 20:49:59 greg Exp $"; +static const char RCSid[] = "$Id: wrapBSDF.c,v 2.16 2015/05/29 07:45:48 greg Exp $"; #endif /* * Wrap BSDF data in valid WINDOW XML file @@ -18,7 +18,6 @@ const char def_template[] = "minimalBSDFt.xml"; const char win6_template[] = "WINDOW6BSDFt.xml"; const char stdin_name[] = ""; - /* input files (can be stdin_name) */ const char *xml_input = NULL; /* unit for materials & geometry */ @@ -26,13 +25,23 @@ const char *attr_unit = "meter"; const char legal_units[] = "meter|foot|inch|centimeter|millimeter"; /* system materials & geometry */ const char *mgf_geometry = NULL; - - /* angle basis */ + /* comment list */ +#define MAXCOMM 30 +const char *commlist[MAXCOMM]; +int ncomm = 0; + /* angle bases */ enum { ABdefault=-1, ABklemsFull=0, ABklemsHalf, ABklemsQuarter, ABtensorTree3, ABtensorTree4, ABend }; int angle_basis = ABdefault; +int correct_solid_angle = 0; + +const char *klems_basis_name[] = { + "LBNL/Klems Full", + "LBNL/Klems Half", + "LBNL/Klems Quarter", +}; /* field IDs and nicknames */ struct s_fieldID { char nickName[4]; @@ -57,11 +66,11 @@ struct s_fieldID { const char *field_assignment[MAXASSIGN]; int nfield_assign = 0; #define FASEP ';' - /* data file(s) & spectra */ enum { DTtransForward, DTtransBackward, DTreflForward, DTreflBackward }; -enum { DSsolar=-1, DSnir=-2, DSxbar31=-3, DSvisible=-4, DSzbar31=-5 }; +enum { DSsolar=-1, DSnir=-2, DSxbar31=-3, DSvisible=-4, DSzbar31=-5, + DSuprime=-6, DSvprime=-7 }; #define MAXFILES 20 @@ -79,24 +88,34 @@ const char top_level_name[] = "WindowElement"; static char basis_definition[][256] = { + "\t\n" "\t\tColumns\n" "\t\t\n" "\t\t\tLBNL/Klems Full\n" - "\t\t\n", + "\t\t\t\n" + "\t\n", + "\t\n" "\t\tColumns\n" "\t\t\n" "\t\t\tLBNL/Klems Half\n" - "\t\t\n", + "\t\t\t\n" + "\t\n", + "\t\n" "\t\tColumns\n" "\t\t\n" "\t\t\tLBNL/Klems Quarter\n" - "\t\t\n", + "\t\t\t\n" + "\t\n", - "\t\tTensorTree3\n", + "\t\n" + "\t\tTensorTree3\n" + "\t\n", - "\t\tTensorTree4\n", + "\t\n" + "\t\tTensorTree4\n" + "\t\n", }; /* Copy data from file descriptor to stdout and close */ @@ -144,6 +163,7 @@ input2str(const char *inpspec) fprintf(stderr, "%s: cannot open\n", inpspec); return ""; } +#ifndef _WIN32 /* XXX somehow broken on Windows */ len = lseek(fd, 0L, SEEK_END); if (len > 0) { lseek(fd, 0L, SEEK_SET); @@ -162,6 +182,7 @@ input2str(const char *inpspec) close(fd); return str; } +#endif fp = fdopen(fd, "r"); /* not a regular file */ } /* reading from stream */ @@ -211,6 +232,7 @@ mat_assignments(const char *caller, const char *fn, ez for (i = 0; i < nfield_assign; i++) { const char *fnext = field_assignment[i]; for ( ; ; ) { + int added = 0; ezxml_t fld; char sbuf[512]; int j; @@ -243,7 +265,9 @@ mat_assignments(const char *caller, const char *fn, ez fprintf(stderr, "%s: warning - adding tag <%s>\n", fn, sbuf); + ezxml_add_txt(wtl, "\t"); fld = ezxml_add_child_d(wtl, sbuf, strlen(wtl->txt)); + ++added; } if (XMLfieldID[j].has_unit) ezxml_set_attr(fld, "unit", attr_unit); @@ -267,6 +291,8 @@ mat_assignments(const char *caller, const char *fn, ez } sbuf[j] = '\0'; ezxml_set_txt_d(fld, sbuf); + if (added) + ezxml_add_txt(wtl, "\n\t"); fnext += (*fnext == FASEP); } } @@ -274,14 +300,12 @@ mat_assignments(const char *caller, const char *fn, ez if (xml_input == win6_template) for (i = 0; XMLfieldID[i].nickName[0]; i++) if (XMLfieldID[i].win_need && - !ezxml_txt(ezxml_child(wtl,XMLfieldID[i].fullName))[0]) { + !ezxml_txt(ezxml_child(wtl,XMLfieldID[i].fullName))[0]) fprintf(stderr, - "%s: missing required '%s' assignment for WINDOW <%s>\n", + "%s: warning - missing '%s' assignment for WINDOW <%s>\n", caller, XMLfieldID[i].nickName, XMLfieldID[i].fullName); - return 0; - } - return 1; /* no errors */ + return 1; } /* Complete angle basis specification */ @@ -309,14 +333,22 @@ finish_angle_basis(ezxml_t ab) sprintf(buf, "%g", i ? .5*(abase_list[n].lat[i].tmin + abase_list[n].lat[i+1].tmin) : .0); + ezxml_add_txt(abb, "\n\t\t\t\t"); ezxml_set_txt_d(ezxml_add_child(abb,"Theta",strlen(abb->txt)), buf); sprintf(buf, "%d", abase_list[n].lat[i].nphis); + ezxml_add_txt(abb, "\n\t\t\t\t"); ezxml_set_txt_d(ezxml_add_child(abb,"nPhis",strlen(abb->txt)), buf); + ezxml_add_txt(abb, "\n\t\t\t\t"); tb = ezxml_add_child(abb, "ThetaBounds", strlen(abb->txt)); + ezxml_add_txt(tb, "\n\t\t\t\t\t"); sprintf(buf, "%g", abase_list[n].lat[i].tmin); - ezxml_set_txt(ezxml_add_child(tb,"LowerTheta",strlen(tb->txt)), buf); + ezxml_set_txt_d(ezxml_add_child(tb,"LowerTheta",strlen(tb->txt)), buf); + ezxml_add_txt(tb, "\n\t\t\t\t\t"); sprintf(buf, "%g", abase_list[n].lat[i+1].tmin); - ezxml_set_txt(ezxml_add_child(tb,"UpperTheta",strlen(tb->txt)), buf); + ezxml_set_txt_d(ezxml_add_child(tb,"UpperTheta",strlen(tb->txt)), buf); + ezxml_add_txt(tb, "\n\t\t\t\t"); + ezxml_add_txt(abb, "\n\t\t\t"); + ezxml_add_txt(ab, "\n\t\t\t"); } return 1; } @@ -332,7 +364,7 @@ determine_angle_basis(const char *fn, ezxml_t wtl) if (wtl == NULL) return -1; ids = ezxml_txt(ezxml_child(wtl, "IncidentDataStructure")); - if (ids == NULL) + if (!ids[0]) return -1; for (i = 0; i < ABend; i++) { ezxml_t parsed = ezxml_parse_str(basis_definition[i], @@ -355,10 +387,55 @@ determine_angle_basis(const char *fn, ezxml_t wtl) return -1; } +/* Filter Klems angle basis, factoring out incident projected solid angle */ +static int +filter_klems_matrix(FILE *fp) +{ +#define MAX_COLUMNS 145 + const char *bn = klems_basis_name[angle_basis]; + float col_corr[MAX_COLUMNS]; + int i, j, n = nabases; + /* get angle basis */ + while (n-- > 0) + if (!strcasecmp(bn, abase_list[n].name)) + break; + if (n < 0) + return 0; + if (abase_list[n].nangles > MAX_COLUMNS) { + fputs("Internal error - too many Klems columns!\n", stderr); + return 0; + } + /* get correction factors */ + for (j = abase_list[n].nangles; j--; ) + col_corr[j] = 1.f / io_getohm(j, &abase_list[n]); + /* read/correct/write matrix */ + for (i = 0; i < abase_list[n].nangles; i++) { + for (j = 0; j < abase_list[n].nangles; j++) { + double d; + if (fscanf(fp, "%lf", &d) != 1) + return 0; + if (d < -1e-3) { + fputs("Negative BSDF data!\n", stderr); + return 0; + } + printf(" %.3e", d*col_corr[j]*(d > 0)); + } + fputc('\n', stdout); + } + while ((i = getc(fp)) != EOF) + if (!isspace(i)) { + fputs("Unexpected data past EOF\n", stderr); + return 0; + } + return 1; /* all is good */ +#undef MAX_COLUMNS +} + /* Write out BSDF data block with surrounding tags */ static int writeBSDFblock(const char *caller, struct s_dfile *df) { + int correct_klems = correct_solid_angle; char *cp; puts("\t"); @@ -366,27 +443,35 @@ writeBSDFblock(const char *caller, struct s_dfile *df) switch (df->spectrum) { case DSvisible: puts("\t\tVisible"); - puts("\t\tSourceSpectrum>CIE Illuminant D65 1nm.ssp"); + puts("\t\tCIE Illuminant D65 1nm.ssp"); puts("\t\tASTM E308 1931 Y.dsp"); break; case DSxbar31: puts("\t\tCIE-X"); - puts("\t\tSourceSpectrum>CIE Illuminant D65 1nm.ssp"); + puts("\t\tCIE Illuminant D65 1nm.ssp"); puts("\t\tASTM E308 1931 X.dsp"); break; case DSzbar31: puts("\t\tCIE-Z"); - puts("\t\tSourceSpectrum>CIE Illuminant D65 1nm.ssp"); + puts("\t\tCIE Illuminant D65 1nm.ssp"); puts("\t\tASTM E308 1931 Z.dsp"); break; + case DSuprime: + puts("\t\tCIE-u"); + puts("\t\tCIE Illuminant D65 1nm.ssp"); + break; + case DSvprime: + puts("\t\tCIE-v"); + puts("\t\tCIE Illuminant D65 1nm.ssp"); + break; case DSsolar: puts("\t\tSolar"); - puts("\t\tSourceSpectrum>CIE Illuminant D65 1nm.ssp"); + puts("\t\tCIE Illuminant D65 1nm.ssp"); puts("\t\tNone"); break; case DSnir: puts("\t\tNIR"); - puts("\t\tSourceSpectrum>PLACE_HOLDER"); + puts("\t\tPLACE_HOLDER"); puts("\t\tPLACE_HOLDER"); break; default: @@ -397,7 +482,7 @@ writeBSDFblock(const char *caller, struct s_dfile *df) spectr_file[df->spectrum]); if (cp != NULL) *cp = '.'; - puts("\t\tSourceSpectrum>CIE Illuminant D65 1nm.ssp"); + puts("\t\tCIE Illuminant D65 1nm.ssp"); printf("\t\t%s\n", spectr_file[df->spectrum]); break; @@ -424,17 +509,19 @@ writeBSDFblock(const char *caller, struct s_dfile *df) puts(""); switch (angle_basis) { case ABklemsFull: - puts("\t\t\tLBNL/Klems Full"); - break; case ABklemsHalf: - puts("\t\t\tLBNL/Klems Half"); - break; case ABklemsQuarter: - puts("\t\t\tLBNL/Klems Quarter"); + fputs("\t\t\t", stdout); + fputs(klems_basis_name[angle_basis], stdout); + puts(""); + fputs("\t\t\t", stdout); + fputs(klems_basis_name[angle_basis], stdout); + puts(""); break; case ABtensorTree3: case ABtensorTree4: puts("\t\t\tLBNL/Shirley-Chiu"); + correct_klems = 0; break; default: fprintf(stderr, "%s: bad angle basis (%d)\n", caller, angle_basis); @@ -443,7 +530,30 @@ writeBSDFblock(const char *caller, struct s_dfile *df) puts("\t\t\tBTDF"); puts("\t\t\t"); fflush(stdout); - if (df->fname == stdin_name) { + if (correct_klems) { /* correct Klems matrix data */ + FILE *fp = stdin; + if (df->fname[0] == '!') + fp = popen(df->fname+1, "r"); + else if (df->fname != stdin_name) + fp = fopen(df->fname, "r"); + if (fp == NULL) { + fprintf(stderr, "%s: cannot open '%s'\n", + caller, df->fname); + return 0; + } + if (!filter_klems_matrix(fp)) { + fprintf(stderr, "%s: Klems data error from '%s'\n", + caller, df->fname); + return 0; + } + if (df->fname[0] != '!') { + fclose(fp); + } else if (pclose(fp)) { + fprintf(stderr, "%s: error running '%s'\n", + caller, df->fname); + return 0; + } + } else if (df->fname == stdin_name) { copy_and_close(fileno(stdin)); } else if (df->fname[0] != '!') { if (!copy_and_close(open(df->fname, O_RDONLY))) { @@ -478,6 +588,9 @@ writeBSDF(const char *caller, ezxml_t fl) free(xml); return 0; } + puts(""); + for (i = 0; i < ncomm; i++) + printf("\n", commlist[i]); fflush(stdout); /* write previous XML info. */ if (write(fileno(stdout), xml, ei) != ei) { free(xml); @@ -490,6 +603,7 @@ writeBSDF(const char *caller, ezxml_t fl) } fputs(xml+ei, stdout); /* write trailer */ free(xml); /* free string */ + fputc('\n', stdout); return (fflush(stdout) == 0); } @@ -570,20 +684,20 @@ wrapBSDF(const char *caller) } /* check basis */ if (angle_basis != ABdefault) { + size_t offset = 0; ezxml_t ab, dd = ezxml_child(wtl, "DataDefinition"); - if (dd == NULL) { - dd = ezxml_add_child(wtl, "DataDefinition", strlen(wtl->txt)); - } else if (dd->child != NULL) { - fprintf(stderr, + if (dd != NULL) { + offset = dd->off; + if (dd->child != NULL) + fprintf(stderr, "%s: warning - replacing existing in '%s'\n", caller, xml_path); - do - ezxml_remove(dd->child); - while (dd->child != NULL); - } - ezxml_insert(ezxml_parse_str(basis_definition[angle_basis], + ezxml_remove(dd); + } else + offset = strlen(wtl->txt); + dd = ezxml_insert(ezxml_parse_str(basis_definition[angle_basis], strlen(basis_definition[angle_basis])), - dd, 0); + wtl, offset); if ((ab = ezxml_child(dd, "AngleBasis")) != NULL && !finish_angle_basis(ab)) goto failure; @@ -608,8 +722,8 @@ UsageExit(const char *pname) { fputs("Usage: ", stderr); fputs(pname, stderr); - fputs(" [-W][-a {kf|kh|kq|t3|t4}][-u unit][-g geom][-f 'x=string;y=string']", stderr); - fputs(" [-s spectr][-tb inp][-tf inp][-rb inp][-rf inp]", stderr); + fputs(" [-W][-c][-a {kf|kh|kq|t3|t4}][-u unit][-g geom][-f 'x=string;y=string']", stderr); + fputs(" [-s spectr][-tb inp][-tf inp][-rb inp][-rf inp][-C comm]", stderr); fputs(" [input.xml]\n", stderr); exit(1); } @@ -675,6 +789,22 @@ main(int argc, char *argv[]) else UsageExit(argv[0]); continue; + case 'c': /* correct solid angle */ + correct_solid_angle = 1; + continue; + case 'C': /* comment */ + if (ncomm >= MAXCOMM) { + fprintf(stderr, "%s: too many comments\n", + argv[0]); + return 1; + } + if (strchr(argv[++i], '>') != NULL) { + fprintf(stderr, "%s: illegal character in comment\n", + argv[0]); + return 1; + } + commlist[ncomm++] = argv[i]; + continue; case 't': /* transmission */ if (i >= argc-1) UsageExit(argv[0]); @@ -716,7 +846,8 @@ main(int argc, char *argv[]) argv[i] = (char *)stdin_name; } data_file[ndataf].fname = argv[i]; - data_file[ndataf++].spectrum = cur_spectrum; + data_file[ndataf].spectrum = cur_spectrum; + ndataf++; continue; case 's': /* spectrum name or input file */ if (++i >= argc) @@ -730,6 +861,10 @@ main(int argc, char *argv[]) cur_spectrum = DSxbar31; else if (!strcasecmp(argv[i], "CIE-Z")) cur_spectrum = DSzbar31; + else if (!strcasecmp(argv[i], "CIE-u")) + cur_spectrum = DSuprime; + else if (!strcasecmp(argv[i], "CIE-v")) + cur_spectrum = DSvprime; else if (!strcasecmp(argv[i], "NIR")) cur_spectrum = DSnir; else { @@ -778,8 +913,6 @@ doneOptions: /* get XML input */ } else { xml_input = argv[i]; } - if ((xml_input == win6_template) & (angle_basis == ABdefault)) - angle_basis = ABklemsFull; /* wrap it! */ return !wrapBSDF(argv[0]); }