--- ray/src/util/dctimestep.c 2009/06/17 20:41:47 2.1 +++ ray/src/util/dctimestep.c 2012/02/21 20:09:25 2.19 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: dctimestep.c,v 2.1 2009/06/17 20:41:47 greg Exp $"; +static const char RCSid[] = "$Id: dctimestep.c,v 2.19 2012/02/21 20:09:25 greg Exp $"; #endif /* * Compute time-step result using Daylight Coefficient method. @@ -10,13 +10,15 @@ static const char RCSid[] = "$Id: dctimestep.c,v 2.1 2 #include #include "standard.h" #include "platform.h" +#include "paths.h" #include "color.h" #include "resolu.h" #include "bsdf.h" +#include "bsdf_m.h" char *progname; /* global argv[0] */ -/* Data types for matrix loading */ +/* Data types for file loading */ enum {DTfromHeader, DTascii, DTfloat, DTdouble, DTrgbe, DTxyze}; /* A color coefficient matrix -- vectors have ncols==1 */ @@ -38,7 +40,7 @@ cm_alloc(int nrows, int ncols) CMATRIX *cm; if ((nrows <= 0) | (ncols <= 0)) - return(NULL); + error(USER, "attempt to create empty matrix"); cm = (CMATRIX *)malloc(sizeof(CMATRIX) + 3*sizeof(COLORV)*(nrows*ncols - 1)); if (cm == NULL) @@ -108,6 +110,8 @@ cm_load(const char *fname, int nrows, int ncols, int d CMATRIX *cm; FILE *fp = stdin; + if (ncols <= 0) + error(USER, "Non-positive number of columns"); if (fname == NULL) fname = ""; else if ((fp = fopen(fname, "r")) == NULL) { @@ -125,7 +129,6 @@ cm_load(const char *fname, int nrows, int ncols, int d break; default: error(USER, "unexpected data type in cm_load()"); - return(NULL); } if (nrows <= 0) { /* don't know length? */ int guessrows = 147; /* usually big enough */ @@ -135,15 +138,21 @@ cm_load(const char *fname, int nrows, int ncols, int d long endpos = ftell(fp); long elemsiz = 3*(dtype==DTfloat ? sizeof(float) : sizeof(double)); - guessrows = (endpos - startpos)/elemsiz; + + if ((endpos - startpos) % (ncols*elemsiz)) { + sprintf(errmsg, + "improper length for binary file '%s'", + fname); + error(USER, errmsg); + } + guessrows = (endpos - startpos)/(ncols*elemsiz); if (fseek(fp, startpos, SEEK_SET) < 0) { sprintf(errmsg, "fseek() error on file '%s'", fname); error(SYSTEM, errmsg); } - if ((endpos - startpos) % elemsiz == 0) - nrows = guessrows; /* confident */ + nrows = guessrows; /* we're confident */ } } cm = cm_alloc(guessrows, ncols); @@ -155,12 +164,12 @@ cm_load(const char *fname, int nrows, int ncols, int d int maxrow = (nrows > 0 ? nrows : 32000); int r, c; for (r = 0; r < maxrow; r++) { - if (r >= cm->nrows) /* need more space? */ + if (r >= cm->nrows) /* need more space? */ cm = cm_resize(cm, 2*cm->nrows); for (c = 0; c < ncols; c++) { COLORV *cv = cm_lval(cm,r,c); if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3) - if ((nrows <= 0) & (r > 0) & (c == 0)) { + if ((nrows <= 0) & (r > 0) & !c) { cm = cm_resize(cm, maxrow=r); break; } else @@ -188,7 +197,7 @@ cm_load(const char *fname, int nrows, int ncols, int d if (nread == cm->nrows*cm->ncols) /* need more space? */ cm = cm_resize(cm, 2*cm->nrows); - else if (nread % cm->ncols == 0) + else if (nread && !(nread % cm->ncols)) /* seem to be done */ cm = cm_resize(cm, nread/cm->ncols); else /* ended mid-row */ @@ -210,7 +219,6 @@ cm_load(const char *fname, int nrows, int ncols, int d copycolor(cvp, dc); cvp += 3; } - } else /* dtype == DTfloat */ { float fc[3]; /* load from float */ COLORV *cvp = cm->cmem; @@ -285,23 +293,79 @@ cm_print(const CMATRIX *cm, FILE *fp) } } -/* convert a BSDF to our matrix representation */ +/* Convert a BSDF to our matrix representation */ static CMATRIX * -cm_bsdf(const struct BSDF_data *bsdf) +cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf) { CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc); - COLORV *mp = cm->cmem; + int nbadohm = 0; + int nneg = 0; int r, c; - - for (r = 0; r < cm->nrows; r++) - for (c = 0; c < cm->ncols; c++, mp += 3) - mp[0] = mp[1] = mp[2] = BSDF_value(bsdf,c,r); + /* reciprocity is "transparent" */ + for (c = 0; c < cm->ncols; c++) { + const double dom = mBSDF_incohm(bsdf,c); + /* projected solid angle */ + nbadohm += (dom <= 0); + + for (r = 0; r < cm->nrows; r++) { + float f = mBSDF_value(bsdf,c,r); + COLORV *mp = cm_lval(cm,r,c); + /* check BSDF value */ + if ((f <= 0) | (dom <= 0)) { + nneg += (f < -FTINY); + f = .0f; + } + copycolor(mp, specCol); + scalecolor(mp, f); + addcolor(mp, bsdfLamb); + scalecolor(mp, dom); + } + } + if (nneg | nbadohm) { + sprintf(errmsg, + "BTDF has %d negatives and %d bad incoming solid angles", + nneg, nbadohm); + error(WARNING, errmsg); + } return(cm); } +/* Load and convert a matrix BSDF from the given XML file */ +static CMATRIX * +cm_loadBSDF(char *fname) +{ + CMATRIX *Tmat; + char *fpath; + SDError ec; + SDData myBSDF; + COLOR bsdfLamb, specCol; + /* find path to BSDF file */ + fpath = getpath(fname, getrlibpath(), R_OK); + if (fpath == NULL) { + sprintf(errmsg, "cannot find BSDF file '%s'", fname); + error(USER, errmsg); + } + SDclearBSDF(&myBSDF, fname); /* load XML and check type */ + ec = SDloadFile(&myBSDF, fpath); + if (ec) + error(USER, transSDError(ec)); + if (myBSDF.tf == NULL || myBSDF.tf->ncomp != 1 || + myBSDF.tf->comp[0].func != &SDhandleMtx) { + sprintf(errmsg, "unsupported BSDF '%s'", fpath); + error(USER, errmsg); + } + /* convert BTDF to matrix */ + ccy2rgb(&myBSDF.tLamb.spec, myBSDF.tLamb.cieY/PI, bsdfLamb); + ccy2rgb(&myBSDF.tf->comp[0].cspec[0], 1., specCol); + Tmat = cm_bsdf(bsdfLamb, specCol, (SDMat *)myBSDF.tf->comp[0].dist); + /* free BSDF and return */ + SDfreeBSDF(&myBSDF); + return(Tmat); +} + /* Sum together a set of images and write result to stdout */ static int -sum_images(const char *fspec, const CMATRIX *cv) +sum_images(const char *fspec, const CMATRIX *cv, FILE *fout) { int myDT = DTfromHeader; CMATRIX *pmat; @@ -317,6 +381,10 @@ sum_images(const char *fspec, const CMATRIX *cv) FILE *fp; int dt, xr, yr; COLORV *psp; + /* check for zero */ + if ((scv[RED] == 0) & (scv[GRN] == 0) & (scv[BLU] == 0) && + (myDT != DTfromHeader) | (i < cv->nrows-1)) + continue; /* open next picture */ sprintf(fname, fspec, i); if ((fp = fopen(fname, "r")) == NULL) { @@ -338,6 +406,11 @@ sum_images(const char *fspec, const CMATRIX *cv) error(SYSTEM, "out of memory in sum_images()"); pmat = cm_alloc(myYR, myXR); memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR); + /* finish header */ + fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, fout); + fputc('\n', fout); + fprtresolu(myXR, myYR, fout); + fflush(fout); } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) { sprintf(errmsg, "picture '%s' format/size mismatch", fname); @@ -360,69 +433,73 @@ sum_images(const char *fspec, const CMATRIX *cv) fclose(fp); /* done this picture */ } free(scanline); - /* finish header */ - fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, stdout); - fputc('\n', stdout); - fprtresolu(myXR, myYR, stdout); /* write scanlines */ for (y = 0; y < myYR; y++) - if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, stdout) < 0) + if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0) return(0); cm_free(pmat); /* all done */ - return(fflush(stdout) == 0); + return(fflush(fout) == 0); } -/* check to see if a string contains a %d specification */ -int -hasDecimalSpec(const char *s) +/* check to see if a string contains a %d or %o specification */ +static int +hasNumberFormat(const char *s) { - while (*s && *s != '%') - s++; - if (!*s) - return(0); - do - ++s; - while (isdigit(*s)); - - return(*s == 'd'); + while (*s) { + while (*s != '%') + if (!*s++) + return(0); + if (*++s == '%') { /* ignore "%%" */ + ++s; + continue; + } + while (isdigit(*s)) /* field length */ + ++s; + /* field we'll use? */ + if ((*s == 'd') | (*s == 'i') | (*s == 'o') | + (*s == 'x') | (*s == 'X')) + return(1); + } + return(0); /* didn't find one */ } int main(int argc, char *argv[]) { - CMATRIX *tvec, *Dmat, *Tmat, *vtmp, *cvec; - struct BSDF_data *btdf; + CMATRIX *cvec; progname = argv[0]; - if ((argc < 4) | (argc > 5)) { - fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n", - progname); + if ((argc < 2) | (argc > 5)) { + fprintf(stderr, "Usage: %s DCspec [tregvec]\n", progname); + fprintf(stderr, " or: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n", + progname); return(1); } - tvec = cm_load(argv[4], 0, 1, DTascii); /* argv[4]==NULL iff argc==4 */ - Dmat = cm_load(argv[3], 0, tvec->nrows, DTfromHeader); - btdf = load_BSDF(argv[2]); - if (btdf == NULL) - return(1); - if (btdf->ninc != Dmat->nrows) { - sprintf(errmsg, "Incoming BTDF dir (%d) mismatch to D (%d)", - btdf->ninc, Dmat->nrows); - error(USER, errmsg); - } + + if (argc > 3) { /* VTDs expression */ + CMATRIX *svec, *Dmat, *Tmat, *ivec; + /* get sky vector */ + svec = cm_load(argv[4], 0, 1, DTascii); + /* load BSDF */ + Tmat = cm_loadBSDF(argv[2]); + /* load Daylight matrix */ + Dmat = cm_load(argv[3], Tmat->ncols, svec->nrows, DTfromHeader); /* multiply vector through */ - vtmp = cm_multiply(Dmat, tvec); - cm_free(Dmat); cm_free(tvec); - Tmat = cm_bsdf(btdf); /* convert BTDF to matrix */ - free_BSDF(btdf); - cvec = cm_multiply(Tmat, vtmp); /* cvec = component vector */ - cm_free(Tmat); cm_free(vtmp); - if (hasDecimalSpec(argv[1])) { /* generating image */ + ivec = cm_multiply(Dmat, svec); + cm_free(Dmat); cm_free(svec); + cvec = cm_multiply(Tmat, ivec); /* cvec = component vector */ + cm_free(Tmat); cm_free(ivec); + } else { /* else just use sky vector */ + cvec = cm_load(argv[2], 0, 1, DTascii); + } + + if (hasNumberFormat(argv[1])) { /* generating image */ SET_FILE_BINARY(stdout); newheader("RADIANCE", stdout); printargs(argc, argv, stdout); fputnow(stdout); - if (!sum_images(argv[1], cvec)) + if (!sum_images(argv[1], cvec, stdout)) return(1); } else { /* generating vector */ CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);