--- ray/src/util/dctimestep.c 2012/08/07 18:17:17 2.20 +++ ray/src/util/dctimestep.c 2013/01/11 17:54:17 2.26 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: dctimestep.c,v 2.20 2012/08/07 18:17:17 greg Exp $"; +static const char RCSid[] = "$Id: dctimestep.c,v 2.26 2013/01/11 17:54:17 greg Exp $"; #endif /* * Compute time-step result using Daylight Coefficient method. @@ -107,8 +107,8 @@ getDTfromHeader(FILE *fp) static CMATRIX * cm_load(const char *fname, int nrows, int ncols, int dtype) { - CMATRIX *cm; FILE *fp = stdin; + CMATRIX *cm; if (ncols <= 0) error(USER, "Non-positive number of columns"); @@ -118,8 +118,11 @@ cm_load(const char *fname, int nrows, int ncols, int d sprintf(errmsg, "cannot open file '%s'", fname); error(SYSTEM, errmsg); } +#ifdef getc_unlocked + flockfile(fp); +#endif if (dtype != DTascii) - SET_FILE_BINARY(fp); + SET_FILE_BINARY(fp); /* doesn't really work */ if (dtype == DTfromHeader) dtype = getDTfromHeader(fp); switch (dtype) { @@ -158,7 +161,7 @@ cm_load(const char *fname, int nrows, int ncols, int d cm = cm_alloc(guessrows, ncols); } else cm = cm_alloc(nrows, ncols); - if (cm == NULL) + if (cm == NULL) /* XXX never happens */ return(NULL); if (dtype == DTascii) { /* read text file */ int maxrow = (nrows > 0 ? nrows : 32000); @@ -233,7 +236,7 @@ cm_load(const char *fname, int nrows, int ncols, int d cvp += 3; } } - if (getc(fp) != EOF) { + if (fgetc(fp) != EOF) { sprintf(errmsg, "unexpected data at end of binary file %s", fname); @@ -242,16 +245,41 @@ cm_load(const char *fname, int nrows, int ncols, int d } if (fp != stdin) fclose(fp); +#ifdef getc_unlocked + else + funlockfile(fp); +#endif return(cm); EOFerror: - sprintf(errmsg, "unexpected EOF reading %s", - fname); + sprintf(errmsg, "unexpected EOF reading %s", fname); error(USER, errmsg); not_handled: error(INTERNAL, "unhandled data size or length in cm_load()"); return(NULL); /* gratis return */ } +/* Extract a column vector from a matrix */ +static CMATRIX * +cm_column(const CMATRIX *cm, int c) +{ + CMATRIX *cvr; + int dr; + + if ((c < 0) | (c >= cm->ncols)) + error(INTERNAL, "column requested outside matrix"); + cvr = cm_alloc(cm->nrows, 1); + if (cvr == NULL) + return(NULL); + for (dr = 0; dr < cm->nrows; dr++) { + const COLORV *sp = cm_lval(cm,dr,c); + COLORV *dp = cv_lval(cvr,dr); + dp[0] = sp[0]; + dp[1] = sp[1]; + dp[2] = sp[2]; + } + return(cvr); +} + /* Scale a matrix by a single value */ static CMATRIX * cm_scale(const CMATRIX *cm1, const COLOR sca) @@ -322,7 +350,7 @@ cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, con int nbadohm = 0; int nneg = 0; int r, c; - /* reciprocity is "transparent" */ + /* loop over incident angles */ for (c = 0; c < cm->ncols; c++) { const double dom = mBSDF_incohm(bsdf,c); /* projected solid angle */ @@ -351,15 +379,80 @@ cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, con return(cm); } +/* Convert between input and output indices for reciprocity */ +static int +recip_out_from_in(const SDMat *bsdf, int in_recip) +{ + FVECT v; + + if (!mBSDF_incvec(v, bsdf, in_recip+.5)) + return(in_recip); /* XXX should be error! */ + v[2] = -v[2]; + return(mBSDF_outndx(bsdf, v)); +} + +/* Convert between output and input indices for reciprocity */ +static int +recip_in_from_out(const SDMat *bsdf, int out_recip) +{ + FVECT v; + + if (!mBSDF_outvec(v, bsdf, out_recip+.5)) + return(out_recip); /* XXX should be error! */ + v[2] = -v[2]; + return(mBSDF_incndx(bsdf, v)); +} + +/* Convert a BSDF to our matrix representation, applying reciprocity */ +static CMATRIX * +cm_bsdf_recip(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf) +{ + CMATRIX *cm = cm_alloc(bsdf->ninc, bsdf->nout); + int nbadohm = 0; + int nneg = 0; + int r, c; + /* loop over incident angles */ + for (c = 0; c < cm->ncols; c++) { + const int ro = recip_out_from_in(bsdf,c); + const double dom = mBSDF_outohm(bsdf,ro); + /* projected solid angle */ + nbadohm += (dom <= 0); + + for (r = 0; r < cm->nrows; r++) { + const int ri = recip_in_from_out(bsdf,r); + float f = mBSDF_value(bsdf,ri,ro); + 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, COLOR cLamb) { - CMATRIX *Tmat; - char *fpath; - SDError ec; - SDData myBSDF; - COLOR bsdfLamb, specCol; + CMATRIX *Tmat; + char *fpath; + int recip; + SDError ec; + SDData myBSDF; + SDSpectralDF *tdf; + COLOR bsdfLamb, specCol; /* find path to BSDF file */ fpath = getpath(fname, getrlibpath(), R_OK); if (fpath == NULL) { @@ -371,19 +464,22 @@ cm_loadBSDF(char *fname, COLOR cLamb) if (ec) error(USER, transSDError(ec)); ccy2rgb(&myBSDF.tLamb.spec, myBSDF.tLamb.cieY/PI, bsdfLamb); - if (myBSDF.tf == NULL) { /* no non-Lambertian transmission? */ + recip = (myBSDF.tb == NULL); + tdf = recip ? myBSDF.tf : myBSDF.tb; + if (tdf == NULL) { /* no non-Lambertian transmission? */ if (cLamb != NULL) copycolor(cLamb, bsdfLamb); SDfreeBSDF(&myBSDF); return(NULL); } - if (myBSDF.tf->ncomp != 1 || myBSDF.tf->comp[0].func != &SDhandleMtx) { + if (tdf->ncomp != 1 || tdf->comp[0].func != &SDhandleMtx) { sprintf(errmsg, "unsupported BSDF '%s'", fpath); error(USER, errmsg); } /* convert BTDF to matrix */ - ccy2rgb(&myBSDF.tf->comp[0].cspec[0], 1., specCol); - Tmat = cm_bsdf(bsdfLamb, specCol, (SDMat *)myBSDF.tf->comp[0].dist); + ccy2rgb(&tdf->comp[0].cspec[0], 1., specCol); + Tmat = recip ? cm_bsdf_recip(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist) + : cm_bsdf(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist); if (cLamb != NULL) /* Lambertian is included */ setcolor(cLamb, .0, .0, .0); /* free BSDF and return */ @@ -391,14 +487,14 @@ cm_loadBSDF(char *fname, COLOR cLamb) return(Tmat); } -/* Sum together a set of images and write result to stdout */ +/* Sum together a set of images and write result to fout */ static int sum_images(const char *fspec, const CMATRIX *cv, FILE *fout) { int myDT = DTfromHeader; - CMATRIX *pmat; - COLOR *scanline; - int myXR, myYR; + COLOR *scanline = NULL; + CMATRIX *pmat = NULL; + int myXR=0, myYR=0; int i, y; if (cv->ncols != 1) @@ -473,6 +569,9 @@ sum_images(const char *fspec, const CMATRIX *cv, FILE static int hasNumberFormat(const char *s) { + if (s == NULL) + return(0); + while (*s) { while (*s != '%') if (!*s++) @@ -494,55 +593,159 @@ hasNumberFormat(const char *s) int main(int argc, char *argv[]) { - CMATRIX *cvec; /* component vector result */ + int skyfmt = DTascii; + int nsteps = 1; + char *ofspec = NULL; + FILE *ofp = stdout; + CMATRIX *cmtx; /* component vector/matrix result */ + char fnbuf[256]; + int a, i; progname = argv[0]; + /* get options */ + for (a = 1; a < argc && argv[a][0] == '-'; a++) + switch (argv[a][1]) { + case 'n': + nsteps = atoi(argv[++a]); + if (nsteps <= 0) + goto userr; + break; + case 'o': + ofspec = argv[++a]; + break; + case 'i': + switch (argv[a][2]) { + case 'f': + skyfmt = DTfloat; + break; + case 'd': + skyfmt = DTdouble; + break; + case 'a': + skyfmt = DTascii; + break; + default: + goto userr; + } + break; + default: + goto userr; + } + if ((argc-a < 1) | (argc-a > 4)) + goto userr; - 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); - } - - if (argc > 3) { /* VTDs expression */ - CMATRIX *svec, *Dmat, *Tmat, *ivec; + if (argc-a > 2) { /* VTDs expression */ + CMATRIX *smtx, *Dmat, *Tmat, *imtx; COLOR tLamb; - /* get sky vector */ - svec = cm_load(argv[4], 0, 1, DTascii); + /* get sky vector/matrix */ + smtx = cm_load(argv[a+3], 0, nsteps, skyfmt); /* load BSDF */ - Tmat = cm_loadBSDF(argv[2], tLamb); + Tmat = cm_loadBSDF(argv[a+1], tLamb); /* load Daylight matrix */ - Dmat = cm_load(argv[3], Tmat==NULL ? 0 : Tmat->ncols, - svec->nrows, DTfromHeader); + Dmat = cm_load(argv[a+2], Tmat==NULL ? 0 : Tmat->ncols, + smtx->nrows, DTfromHeader); /* multiply vector through */ - ivec = cm_multiply(Dmat, svec); - cm_free(Dmat); cm_free(svec); + imtx = cm_multiply(Dmat, smtx); + cm_free(Dmat); cm_free(smtx); if (Tmat == NULL) { /* diffuse only */ - cvec = cm_scale(ivec, tLamb); + cmtx = cm_scale(imtx, tLamb); } else { /* else apply BTDF matrix */ - cvec = cm_multiply(Tmat, ivec); + cmtx = cm_multiply(Tmat, imtx); cm_free(Tmat); } - cm_free(ivec); - } else { /* else just use sky vector */ - cvec = cm_load(argv[2], 0, 1, DTascii); + cm_free(imtx); + } else { /* sky vector/matrix only */ + cmtx = cm_load(argv[a+1], 0, nsteps, skyfmt); } - - 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, stdout)) + /* prepare output stream */ + if ((ofspec != NULL) & (nsteps == 1) && hasNumberFormat(ofspec)) { + sprintf(fnbuf, ofspec, 1); + ofspec = fnbuf; + } + if (ofspec != NULL && !hasNumberFormat(ofspec)) { + if ((ofp = fopen(ofspec, "w")) == NULL) { + fprintf(stderr, "%s: cannot open '%s' for output\n", + progname, ofspec); return(1); - } else { /* generating vector */ - CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader); - CMATRIX *rvec = cm_multiply(Vmat, cvec); + } + ofspec = NULL; /* only need to open once */ + } + if (hasNumberFormat(argv[a])) { /* generating image(s) */ + if (ofspec == NULL) { + SET_FILE_BINARY(ofp); + newheader("RADIANCE", ofp); + printargs(argc, argv, ofp); + fputnow(ofp); + } + if (nsteps > 1) /* multiple output frames? */ + for (i = 0; i < nsteps; i++) { + CMATRIX *cvec = cm_column(cmtx, i); + if (ofspec != NULL) { + sprintf(fnbuf, ofspec, i+1); + if ((ofp = fopen(fnbuf, "wb")) == NULL) { + fprintf(stderr, + "%s: cannot open '%s' for output\n", + progname, fnbuf); + return(1); + } + newheader("RADIANCE", ofp); + printargs(argc, argv, ofp); + fputnow(ofp); + } + fprintf(ofp, "FRAME=%d\n", i+1); + if (!sum_images(argv[a], cvec, ofp)) + return(1); + if (ofspec != NULL) { + if (fclose(ofp) == EOF) { + fprintf(stderr, + "%s: error writing to '%s'\n", + progname, fnbuf); + return(1); + } + ofp = stdout; + } + cm_free(cvec); + } + else if (!sum_images(argv[a], cmtx, ofp)) + return(1); + } else { /* generating vector/matrix */ + CMATRIX *Vmat = cm_load(argv[a], 0, cmtx->nrows, DTfromHeader); + CMATRIX *rmtx = cm_multiply(Vmat, cmtx); cm_free(Vmat); - cm_print(rvec, stdout); - cm_free(rvec); + if (ofspec != NULL) /* multiple vector files? */ + for (i = 0; i < nsteps; i++) { + CMATRIX *rvec = cm_column(rmtx, i); + sprintf(fnbuf, ofspec, i+1); + if ((ofp = fopen(fnbuf, "w")) == NULL) { + fprintf(stderr, + "%s: cannot open '%s' for output\n", + progname, fnbuf); + return(1); + } + cm_print(rvec, ofp); + if (fclose(ofp) == EOF) { + fprintf(stderr, + "%s: error writing to '%s'\n", + progname, fnbuf); + return(1); + } + ofp = stdout; + cm_free(rvec); + } + else + cm_print(rmtx, ofp); + cm_free(rmtx); } - cm_free(cvec); /* final clean-up */ + if (fflush(ofp) == EOF) { /* final clean-up */ + fprintf(stderr, "%s: write error on output\n", progname); + return(1); + } + cm_free(cmtx); return(0); +userr: + fprintf(stderr, "Usage: %s [-n nsteps][-o ospec][-i{f|d}] DCspec [skyf]\n", + progname); + fprintf(stderr, " or: %s [-n nsteps][-o ospec][-i{f|d}] Vspec Tbsdf.xml Dmat.dat [skyf]\n", + progname); + return(1); }