src/util/dctimestep.c

#ifndef lint
static const char RCSid[] = "$Id: dctimestep.c,v 2.25 2013/01/11 17:21:39 greg Exp $";
#endif
/*
 * Compute time-step result using Daylight Coefficient method.
 *
 *      G. Ward
 */

#include <ctype.h>
#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 file loading */
enum {DTfromHeader, DTascii, DTfloat, DTdouble, DTrgbe, DTxyze};

/* A color coefficient matrix -- vectors have ncols==1 */
typedef struct {
        int     nrows, ncols;
        COLORV  cmem[3];                /* extends struct */
} CMATRIX;

#define COLSPEC (sizeof(COLORV)==sizeof(float) ? "%f %f %f" : "%lf %lf %lf")

#define cm_lval(cm,r,c) ((cm)->cmem + 3*((r)*(cm)->ncols + (c)))

#define cv_lval(cm,i)   ((cm)->cmem + 3*(i))

/* Allocate a color coefficient matrix */
static CMATRIX *
cm_alloc(int nrows, int ncols)
{
        CMATRIX *cm;

        if ((nrows <= 0) | (ncols <= 0))
                error(USER, "attempt to create empty matrix");
        cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
                                3*sizeof(COLORV)*(nrows*ncols - 1));
        if (cm == NULL)
                error(SYSTEM, "out of memory in cm_alloc()");
        cm->nrows = nrows;
        cm->ncols = ncols;
        return(cm);
}

#define cm_free(cm)     free(cm)

/* Resize color coefficient matrix */
static CMATRIX *
cm_resize(CMATRIX *cm, int nrows)
{
        if (nrows == cm->nrows)
                return(cm);
        if (nrows <= 0) {
                cm_free(cm);
                return(NULL);
        }
        cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) +
                        3*sizeof(COLORV)*(nrows*cm->ncols - 1));
        if (cm == NULL)
                error(SYSTEM, "out of memory in cm_resize()");
        cm->nrows = nrows;
        return(cm);
}

/* Load header to obtain data type */
static int
getDT(char *s, void *p)
{
        char    fmt[32];
        
        if (formatval(fmt, s)) {
                if (!strcmp(fmt, "ascii"))
                        *((int *)p) = DTascii;
                else if (!strcmp(fmt, "float"))
                        *((int *)p) = DTfloat;
                else if (!strcmp(fmt, "double"))
                        *((int *)p) = DTdouble;
                else if (!strcmp(fmt, COLRFMT))
                        *((int *)p) = DTrgbe;
                else if (!strcmp(fmt, CIEFMT))
                        *((int *)p) = DTxyze;
        }
        return(0);
}

static int
getDTfromHeader(FILE *fp)
{
        int     dt = DTfromHeader;
        
        if (getheader(fp, getDT, &dt) < 0)
                error(SYSTEM, "header read error");
        if (dt == DTfromHeader)
                error(USER, "missing data format in header");
        return(dt);
}

/* Allocate and load a matrix from the given file (or stdin if NULL) */
static CMATRIX *
cm_load(const char *fname, int nrows, int ncols, int dtype)
{
        FILE    *fp = stdin;
        CMATRIX *cm;

        if (ncols <= 0)
                error(USER, "Non-positive number of columns");
        if (fname == NULL)
                fname = "<stdin>";
        else if ((fp = fopen(fname, "r")) == NULL) {
                sprintf(errmsg, "cannot open file '%s'", fname);
                error(SYSTEM, errmsg);
        }
#ifdef getc_unlocked
        flockfile(fp);
#endif
        if (dtype != DTascii)
                SET_FILE_BINARY(fp);            /* doesn't really work */
        if (dtype == DTfromHeader)
                dtype = getDTfromHeader(fp);
        switch (dtype) {
        case DTascii:
        case DTfloat:
        case DTdouble:
                break;
        default:
                error(USER, "unexpected data type in cm_load()");
        }
        if (nrows <= 0) {                       /* don't know length? */
                int     guessrows = 147;        /* usually big enough */
                if ((dtype != DTascii) & (fp != stdin)) {
                        long    startpos = ftell(fp);
                        if (fseek(fp, 0L, SEEK_END) == 0) {
                                long    endpos = ftell(fp);
                                long    elemsiz = 3*(dtype==DTfloat ?
                                            sizeof(float) : sizeof(double));

                                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);
                                }
                                nrows = guessrows;      /* we're confident */
                        }
                }
                cm = cm_alloc(guessrows, ncols);
        } else
                cm = cm_alloc(nrows, ncols);
        if (cm == NULL)                                 /* XXX never happens */
                return(NULL);
        if (dtype == DTascii) {                         /* read text file */
                int     maxrow = (nrows > 0 ? nrows : 32000);
                int     r, c;
                for (r = 0; r < maxrow; r++) {
                    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) {
                                        cm = cm_resize(cm, maxrow=r);
                                        break;
                                } else
                                        goto EOFerror;
                    }
                }
                while ((c = getc(fp)) != EOF)
                        if (!isspace(c)) {
                                sprintf(errmsg,
                                "unexpected data at end of ascii file %s",
                                                fname);
                                error(WARNING, errmsg);
                                break;
                        }
        } else {                                        /* read binary file */
                if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) :
                                                        sizeof(double))) {
                        int     nread = 0;
                        do {                            /* read all we can */
                                nread += fread(cm->cmem + 3*nread,
                                                3*sizeof(COLORV),
                                                cm->nrows*cm->ncols - nread,
                                                fp);
                                if (nrows <= 0) {       /* unknown length */
                                        if (nread == cm->nrows*cm->ncols)
                                                        /* need more space? */
                                                cm = cm_resize(cm, 2*cm->nrows);
                                        else if (nread && !(nread % cm->ncols))
                                                        /* seem to be  done */
                                                cm = cm_resize(cm, nread/cm->ncols);
                                        else            /* ended mid-row */
                                                goto EOFerror;
                                } else if (nread < cm->nrows*cm->ncols)
                                        goto EOFerror;
                        } while (nread < cm->nrows*cm->ncols);

                } else if (dtype == DTdouble) {
                        double  dc[3];                  /* load from double */
                        COLORV  *cvp = cm->cmem;
                        int     n = nrows*ncols;

                        if (n <= 0)
                                goto not_handled;
                        while (n--) {
                                if (fread(dc, sizeof(double), 3, fp) != 3)
                                        goto EOFerror;
                                copycolor(cvp, dc);
                                cvp += 3;
                        }
                } else /* dtype == DTfloat */ {
                        float   fc[3];                  /* load from float */
                        COLORV  *cvp = cm->cmem;
                        int     n = nrows*ncols;

                        if (n <= 0)
                                goto not_handled;
                        while (n--) {
                                if (fread(fc, sizeof(float), 3, fp) != 3)
                                        goto EOFerror;
                                copycolor(cvp, fc);
                                cvp += 3;
                        }
                }
                if (fgetc(fp) != EOF) {
                                sprintf(errmsg,
                                "unexpected data at end of binary file %s",
                                                fname);
                                error(WARNING, errmsg);
                }
        }
        if (fp != stdin)
                fclose(fp);
#ifdef getc_unlocked
        else
                funlockfile(fp);
#endif
        return(cm);
EOFerror:
        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)
{
        CMATRIX *cmr;
        int     dr, dc;

        cmr = cm_alloc(cm1->nrows, cm1->ncols);
        if (cmr == NULL)
                return(NULL);
        for (dr = 0; dr < cmr->nrows; dr++)
            for (dc = 0; dc < cmr->ncols; dc++) {
                const COLORV    *sp = cm_lval(cm1,dr,dc);
                COLORV          *dp = cm_lval(cmr,dr,dc);
                dp[0] = sp[0] * sca[0];
                dp[1] = sp[1] * sca[1];
                dp[2] = sp[2] * sca[2];
            }
        return(cmr);
}

/* Multiply two matrices (or a matrix and a vector) and allocate the result */
static CMATRIX *
cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
{
        CMATRIX *cmr;
        int     dr, dc, i;

        if ((cm1->ncols <= 0) | (cm1->ncols != cm2->nrows))
                error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
        cmr = cm_alloc(cm1->nrows, cm2->ncols);
        if (cmr == NULL)
                return(NULL);
        for (dr = 0; dr < cmr->nrows; dr++)
            for (dc = 0; dc < cmr->ncols; dc++) {
                COLORV  *dp = cm_lval(cmr,dr,dc);
                dp[0] = dp[1] = dp[2] = 0;
                for (i = 0; i < cm1->ncols; i++) {
                    const COLORV        *cp1 = cm_lval(cm1,dr,i);
                    const COLORV        *cp2 = cm_lval(cm2,i,dc);
                    dp[0] += cp1[0] * cp2[0];
                    dp[1] += cp1[1] * cp2[1];
                    dp[2] += cp1[2] * cp2[2];
                }
            }
        return(cmr);
}

/* print out matrix as ASCII text -- no header */
static void
cm_print(const CMATRIX *cm, FILE *fp)
{
        int             r, c;
        const COLORV    *mp = cm->cmem;
        
        for (r = 0; r < cm->nrows; r++) {
                for (c = 0; c < cm->ncols; c++, mp += 3)
                        fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]);
                fputc('\n', fp);
        }
}

/* Convert a BSDF to our matrix representation */
static CMATRIX *
cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf)
{
        CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
        int     nbadohm = 0;
        int     nneg = 0;
        int     r, c;
                                        /* loop over incident angles */
        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);
}

/* 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;
        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) {
                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));
        ccy2rgb(&myBSDF.tLamb.spec, myBSDF.tLamb.cieY/PI, bsdfLamb);
        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 (tdf->ncomp != 1 || tdf->comp[0].func != &SDhandleMtx) {
                sprintf(errmsg, "unsupported BSDF '%s'", fpath);
                error(USER, errmsg);
        }
                                        /* convert BTDF to matrix */
        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 */
        SDfreeBSDF(&myBSDF);
        return(Tmat);
}

/* 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;
        COLOR   *scanline = NULL;
        CMATRIX *pmat = NULL;
        int     myXR=0, myYR=0;
        int     i, y;

        if (cv->ncols != 1)
                error(INTERNAL, "expected vector in sum_images()");
        for (i = 0; i < cv->nrows; i++) {
                const COLORV    *scv = cv_lval(cv,i);
                char            fname[1024];
                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) {
                        sprintf(errmsg, "cannot open picture '%s'", fname);
                        error(SYSTEM, errmsg);
                }
                SET_FILE_BINARY(fp);
                dt = getDTfromHeader(fp);
                if ((dt != DTrgbe) & (dt != DTxyze) ||
                                !fscnresolu(&xr, &yr, fp)) {
                        sprintf(errmsg, "file '%s' not a picture", fname);
                        error(USER, errmsg);
                }
                if (myDT == DTfromHeader) {             /* on first one */
                        myDT = dt;
                        myXR = xr; myYR = yr;
                        scanline = (COLOR *)malloc(sizeof(COLOR)*myXR);
                        if (scanline == NULL)
                                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);
                        error(USER, errmsg);
                }
                psp = pmat->cmem;
                for (y = 0; y < yr; y++) {              /* read it in */
                        int     x;
                        if (freadscan(scanline, xr, fp) < 0) {
                                sprintf(errmsg, "error reading picture '%s'",
                                                fname);
                                error(SYSTEM, errmsg);
                        }
                                                        /* sum in scanline */
                        for (x = 0; x < xr; x++, psp += 3) {
                                multcolor(scanline[x], scv);
                                addcolor(psp, scanline[x]);
                        }
                }
                fclose(fp);                             /* done this picture */
        }
        free(scanline);
                                                        /* write scanlines */
        for (y = 0; y < myYR; y++)
                if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0)
                        return(0);
        cm_free(pmat);                                  /* all done */
        return(fflush(fout) == 0);
}

/* check to see if a string contains a %d or %o specification */
static int
hasNumberFormat(const char *s)
{
        if (s == NULL)
                return(0);

        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[])
{
        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-a > 2) {                       /* VTDs expression */
                CMATRIX *smtx, *Dmat, *Tmat, *imtx;
                COLOR   tLamb;
                                                /* get sky vector/matrix */
                smtx = cm_load(argv[a+3], 0, nsteps, skyfmt);
                                                /* load BSDF */
                Tmat = cm_loadBSDF(argv[a+1], tLamb);
                                                /* load Daylight matrix */
                Dmat = cm_load(argv[a+2], Tmat==NULL ? 0 : Tmat->ncols,
                                        smtx->nrows, DTfromHeader);
                                                /* multiply vector through */
                imtx = cm_multiply(Dmat, smtx);
                cm_free(Dmat); cm_free(smtx);
                if (Tmat == NULL) {             /* diffuse only */
                        cmtx = cm_scale(imtx, tLamb);
                } else {                        /* else apply BTDF matrix */
                        cmtx = cm_multiply(Tmat, imtx);
                        cm_free(Tmat); 
                }
                cm_free(imtx);
        } else {                                /* sky vector/matrix only */
                cmtx = cm_load(argv[a+1], 0, nsteps, skyfmt);
        }
                                                /* 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);
                }
                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);
                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);
        }
        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);
}
Revision:	2.26
Committed:	Fri Jan 11 17:54:17 2013 UTC (12 years, 9 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.25:	+3 -3 lines
Log Message:	Minor fix in unreachable error
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.26	static const char RCSid[] = "$Id: dctimestep.c,v 2.25 2013/01/11 17:21:39 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Compute time-step result using Daylight Coefficient method.
6			*
7			* G. Ward
8			*/
9
10			#include <ctype.h>
11			#include "standard.h"
12			#include "platform.h"
13	greg	2.15	#include "paths.h"
14	greg	2.1	#include "color.h"
15			#include "resolu.h"
16			#include "bsdf.h"
17	greg	2.15	#include "bsdf_m.h"
18	greg	2.1
19			char progname; / global argv[0] */
20
21	greg	2.2	/* Data types for file loading */
22	greg	2.1	enum {DTfromHeader, DTascii, DTfloat, DTdouble, DTrgbe, DTxyze};
23
24			/* A color coefficient matrix -- vectors have ncols==1 */
25			typedef struct {
26			int nrows, ncols;
27			COLORV cmem[3]; /* extends struct */
28			} CMATRIX;
29
30			#define COLSPEC (sizeof(COLORV)==sizeof(float) ? "%f %f %f" : "%lf %lf %lf")
31
32			#define cm_lval(cm,r,c) ((cm)->cmem + 3((r)(cm)->ncols + (c)))
33
34			#define cv_lval(cm,i) ((cm)->cmem + 3*(i))
35
36			/* Allocate a color coefficient matrix */
37			static CMATRIX *
38			cm_alloc(int nrows, int ncols)
39			{
40			CMATRIX *cm;
41
42			if ((nrows <= 0) \| (ncols <= 0))
43	greg	2.10	error(USER, "attempt to create empty matrix");
44	greg	2.1	cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
45			3sizeof(COLORV)(nrows*ncols - 1));
46			if (cm == NULL)
47			error(SYSTEM, "out of memory in cm_alloc()");
48			cm->nrows = nrows;
49			cm->ncols = ncols;
50			return(cm);
51			}
52
53			#define cm_free(cm) free(cm)
54
55			/* Resize color coefficient matrix */
56			static CMATRIX *
57			cm_resize(CMATRIX *cm, int nrows)
58			{
59			if (nrows == cm->nrows)
60			return(cm);
61			if (nrows <= 0) {
62			cm_free(cm);
63			return(NULL);
64			}
65			cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) +
66			3sizeof(COLORV)(nrows*cm->ncols - 1));
67			if (cm == NULL)
68			error(SYSTEM, "out of memory in cm_resize()");
69			cm->nrows = nrows;
70			return(cm);
71			}
72
73			/* Load header to obtain data type */
74			static int
75			getDT(char s, void p)
76			{
77			char fmt[32];
78
79			if (formatval(fmt, s)) {
80			if (!strcmp(fmt, "ascii"))
81			((int )p) = DTascii;
82			else if (!strcmp(fmt, "float"))
83			((int )p) = DTfloat;
84			else if (!strcmp(fmt, "double"))
85			((int )p) = DTdouble;
86			else if (!strcmp(fmt, COLRFMT))
87			((int )p) = DTrgbe;
88			else if (!strcmp(fmt, CIEFMT))
89			((int )p) = DTxyze;
90			}
91			return(0);
92			}
93
94			static int
95			getDTfromHeader(FILE *fp)
96			{
97			int dt = DTfromHeader;
98
99			if (getheader(fp, getDT, &dt) < 0)
100			error(SYSTEM, "header read error");
101			if (dt == DTfromHeader)
102			error(USER, "missing data format in header");
103			return(dt);
104			}
105
106			/* Allocate and load a matrix from the given file (or stdin if NULL) */
107			static CMATRIX *
108			cm_load(const char *fname, int nrows, int ncols, int dtype)
109			{
110	greg	2.25	FILE *fp = stdin;
111	greg	2.1	CMATRIX *cm;
112
113	greg	2.12	if (ncols <= 0)
114			error(USER, "Non-positive number of columns");
115	greg	2.1	if (fname == NULL)
116			fname = "<stdin>";
117			else if ((fp = fopen(fname, "r")) == NULL) {
118			sprintf(errmsg, "cannot open file '%s'", fname);
119			error(SYSTEM, errmsg);
120			}
121	greg	2.25	#ifdef getc_unlocked
122			flockfile(fp);
123			#endif
124	greg	2.1	if (dtype != DTascii)
125	greg	2.26	SET_FILE_BINARY(fp); /* doesn't really work */
126	greg	2.1	if (dtype == DTfromHeader)
127			dtype = getDTfromHeader(fp);
128			switch (dtype) {
129			case DTascii:
130			case DTfloat:
131			case DTdouble:
132			break;
133			default:
134			error(USER, "unexpected data type in cm_load()");
135			}
136			if (nrows <= 0) { /* don't know length? */
137			int guessrows = 147; /* usually big enough */
138			if ((dtype != DTascii) & (fp != stdin)) {
139			long startpos = ftell(fp);
140			if (fseek(fp, 0L, SEEK_END) == 0) {
141			long endpos = ftell(fp);
142			long elemsiz = 3*(dtype==DTfloat ?
143			sizeof(float) : sizeof(double));
144	greg	2.2
145			if ((endpos - startpos) % (ncols*elemsiz)) {
146			sprintf(errmsg,
147			"improper length for binary file '%s'",
148			fname);
149			error(USER, errmsg);
150			}
151			guessrows = (endpos - startpos)/(ncols*elemsiz);
152	greg	2.1	if (fseek(fp, startpos, SEEK_SET) < 0) {
153			sprintf(errmsg,
154			"fseek() error on file '%s'",
155			fname);
156			error(SYSTEM, errmsg);
157			}
158	greg	2.2	nrows = guessrows; /* we're confident */
159	greg	2.1	}
160			}
161			cm = cm_alloc(guessrows, ncols);
162			} else
163			cm = cm_alloc(nrows, ncols);
164	greg	2.25	if (cm == NULL) /* XXX never happens */
165	greg	2.1	return(NULL);
166			if (dtype == DTascii) { /* read text file */
167			int maxrow = (nrows > 0 ? nrows : 32000);
168			int r, c;
169			for (r = 0; r < maxrow; r++) {
170	greg	2.6	if (r >= cm->nrows) /* need more space? */
171	greg	2.1	cm = cm_resize(cm, 2*cm->nrows);
172			for (c = 0; c < ncols; c++) {
173			COLORV *cv = cm_lval(cm,r,c);
174			if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3)
175	greg	2.6	if ((nrows <= 0) & (r > 0) & !c) {
176	greg	2.1	cm = cm_resize(cm, maxrow=r);
177			break;
178			} else
179			goto EOFerror;
180			}
181			}
182			while ((c = getc(fp)) != EOF)
183			if (!isspace(c)) {
184			sprintf(errmsg,
185			"unexpected data at end of ascii file %s",
186			fname);
187			error(WARNING, errmsg);
188			break;
189			}
190			} else { /* read binary file */
191			if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) :
192			sizeof(double))) {
193			int nread = 0;
194			do { /* read all we can */
195			nread += fread(cm->cmem + 3*nread,
196			3*sizeof(COLORV),
197			cm->nrows*cm->ncols - nread,
198			fp);
199			if (nrows <= 0) { /* unknown length */
200			if (nread == cm->nrows*cm->ncols)
201			/* need more space? */
202			cm = cm_resize(cm, 2*cm->nrows);
203	greg	2.6	else if (nread && !(nread % cm->ncols))
204	greg	2.1	/* seem to be done */
205			cm = cm_resize(cm, nread/cm->ncols);
206			else /* ended mid-row */
207			goto EOFerror;
208			} else if (nread < cm->nrows*cm->ncols)
209			goto EOFerror;
210			} while (nread < cm->nrows*cm->ncols);
211
212			} else if (dtype == DTdouble) {
213			double dc[3]; /* load from double */
214			COLORV *cvp = cm->cmem;
215			int n = nrows*ncols;
216
217			if (n <= 0)
218			goto not_handled;
219			while (n--) {
220			if (fread(dc, sizeof(double), 3, fp) != 3)
221			goto EOFerror;
222			copycolor(cvp, dc);
223			cvp += 3;
224			}
225			} else /* dtype == DTfloat */ {
226			float fc[3]; /* load from float */
227			COLORV *cvp = cm->cmem;
228			int n = nrows*ncols;
229
230			if (n <= 0)
231			goto not_handled;
232			while (n--) {
233			if (fread(fc, sizeof(float), 3, fp) != 3)
234			goto EOFerror;
235			copycolor(cvp, fc);
236			cvp += 3;
237			}
238			}
239	greg	2.25	if (fgetc(fp) != EOF) {
240	greg	2.1	sprintf(errmsg,
241			"unexpected data at end of binary file %s",
242			fname);
243			error(WARNING, errmsg);
244			}
245			}
246			if (fp != stdin)
247			fclose(fp);
248	greg	2.25	#ifdef getc_unlocked
249			else
250			funlockfile(fp);
251			#endif
252	greg	2.1	return(cm);
253			EOFerror:
254	greg	2.25	sprintf(errmsg, "unexpected EOF reading %s", fname);
255	greg	2.1	error(USER, errmsg);
256			not_handled:
257			error(INTERNAL, "unhandled data size or length in cm_load()");
258			return(NULL); /* gratis return */
259			}
260
261	greg	2.24	/* Extract a column vector from a matrix */
262			static CMATRIX *
263			cm_column(const CMATRIX *cm, int c)
264			{
265			CMATRIX *cvr;
266			int dr;
267
268			if ((c < 0) \| (c >= cm->ncols))
269	greg	2.26	error(INTERNAL, "column requested outside matrix");
270	greg	2.24	cvr = cm_alloc(cm->nrows, 1);
271			if (cvr == NULL)
272			return(NULL);
273			for (dr = 0; dr < cm->nrows; dr++) {
274			const COLORV *sp = cm_lval(cm,dr,c);
275			COLORV *dp = cv_lval(cvr,dr);
276			dp[0] = sp[0];
277			dp[1] = sp[1];
278			dp[2] = sp[2];
279			}
280			return(cvr);
281			}
282
283	greg	2.20	/* Scale a matrix by a single value */
284			static CMATRIX *
285			cm_scale(const CMATRIX *cm1, const COLOR sca)
286			{
287			CMATRIX *cmr;
288			int dr, dc;
289
290			cmr = cm_alloc(cm1->nrows, cm1->ncols);
291			if (cmr == NULL)
292			return(NULL);
293			for (dr = 0; dr < cmr->nrows; dr++)
294			for (dc = 0; dc < cmr->ncols; dc++) {
295			const COLORV *sp = cm_lval(cm1,dr,dc);
296			COLORV *dp = cm_lval(cmr,dr,dc);
297			dp[0] = sp[0] * sca[0];
298			dp[1] = sp[1] * sca[1];
299			dp[2] = sp[2] * sca[2];
300			}
301			return(cmr);
302			}
303
304			/* Multiply two matrices (or a matrix and a vector) and allocate the result */
305	greg	2.1	static CMATRIX *
306			cm_multiply(const CMATRIX cm1, const CMATRIX cm2)
307			{
308			CMATRIX *cmr;
309			int dr, dc, i;
310
311			if ((cm1->ncols <= 0) \| (cm1->ncols != cm2->nrows))
312			error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
313			cmr = cm_alloc(cm1->nrows, cm2->ncols);
314			if (cmr == NULL)
315			return(NULL);
316			for (dr = 0; dr < cmr->nrows; dr++)
317			for (dc = 0; dc < cmr->ncols; dc++) {
318			COLORV *dp = cm_lval(cmr,dr,dc);
319			dp[0] = dp[1] = dp[2] = 0;
320			for (i = 0; i < cm1->ncols; i++) {
321			const COLORV *cp1 = cm_lval(cm1,dr,i);
322			const COLORV *cp2 = cm_lval(cm2,i,dc);
323			dp[0] += cp1[0] * cp2[0];
324			dp[1] += cp1[1] * cp2[1];
325			dp[2] += cp1[2] * cp2[2];
326			}
327			}
328			return(cmr);
329			}
330
331			/* print out matrix as ASCII text -- no header */
332			static void
333			cm_print(const CMATRIX cm, FILE fp)
334			{
335			int r, c;
336			const COLORV *mp = cm->cmem;
337
338			for (r = 0; r < cm->nrows; r++) {
339			for (c = 0; c < cm->ncols; c++, mp += 3)
340			fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]);
341			fputc('\n', fp);
342			}
343			}
344
345	greg	2.15	/* Convert a BSDF to our matrix representation */
346	greg	2.1	static CMATRIX *
347	greg	2.15	cm_bsdf(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf)
348	greg	2.1	{
349	greg	2.16	CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
350	greg	2.3	int nbadohm = 0;
351			int nneg = 0;
352	greg	2.1	int r, c;
353	greg	2.22	/* loop over incident angles */
354	greg	2.16	for (c = 0; c < cm->ncols; c++) {
355	greg	2.17	const double dom = mBSDF_incohm(bsdf,c);
356			/* projected solid angle */
357	greg	2.16	nbadohm += (dom <= 0);
358
359			for (r = 0; r < cm->nrows; r++) {
360			float f = mBSDF_value(bsdf,c,r);
361			COLORV *mp = cm_lval(cm,r,c);
362			/* check BSDF value */
363			if ((f <= 0) \| (dom <= 0)) {
364			nneg += (f < -FTINY);
365			f = .0f;
366	greg	2.3	}
367	greg	2.16	copycolor(mp, specCol);
368			scalecolor(mp, f);
369			addcolor(mp, bsdfLamb);
370			scalecolor(mp, dom);
371	greg	2.3	}
372	greg	2.7	}
373	greg	2.15	if (nneg \| nbadohm) {
374	greg	2.3	sprintf(errmsg,
375	greg	2.4	"BTDF has %d negatives and %d bad incoming solid angles",
376			nneg, nbadohm);
377	greg	2.3	error(WARNING, errmsg);
378			}
379	greg	2.1	return(cm);
380			}
381
382	greg	2.22	/* Convert between input and output indices for reciprocity */
383			static int
384			recip_out_from_in(const SDMat *bsdf, int in_recip)
385			{
386			FVECT v;
387
388			if (!mBSDF_incvec(v, bsdf, in_recip+.5))
389			return(in_recip); /* XXX should be error! */
390			v[2] = -v[2];
391			return(mBSDF_outndx(bsdf, v));
392			}
393
394			/* Convert between output and input indices for reciprocity */
395			static int
396			recip_in_from_out(const SDMat *bsdf, int out_recip)
397			{
398			FVECT v;
399
400			if (!mBSDF_outvec(v, bsdf, out_recip+.5))
401			return(out_recip); /* XXX should be error! */
402			v[2] = -v[2];
403			return(mBSDF_incndx(bsdf, v));
404			}
405
406			/* Convert a BSDF to our matrix representation, applying reciprocity */
407			static CMATRIX *
408			cm_bsdf_recip(const COLOR bsdfLamb, const COLOR specCol, const SDMat *bsdf)
409			{
410			CMATRIX *cm = cm_alloc(bsdf->ninc, bsdf->nout);
411			int nbadohm = 0;
412			int nneg = 0;
413			int r, c;
414			/* loop over incident angles */
415			for (c = 0; c < cm->ncols; c++) {
416			const int ro = recip_out_from_in(bsdf,c);
417			const double dom = mBSDF_outohm(bsdf,ro);
418			/* projected solid angle */
419			nbadohm += (dom <= 0);
420
421			for (r = 0; r < cm->nrows; r++) {
422			const int ri = recip_in_from_out(bsdf,r);
423			float f = mBSDF_value(bsdf,ri,ro);
424			COLORV *mp = cm_lval(cm,r,c);
425			/* check BSDF value */
426			if ((f <= 0) \| (dom <= 0)) {
427			nneg += (f < -FTINY);
428			f = .0f;
429			}
430			copycolor(mp, specCol);
431			scalecolor(mp, f);
432			addcolor(mp, bsdfLamb);
433			scalecolor(mp, dom);
434			}
435			}
436			if (nneg \| nbadohm) {
437			sprintf(errmsg,
438			"BTDF has %d negatives and %d bad incoming solid angles",
439			nneg, nbadohm);
440			error(WARNING, errmsg);
441			}
442			return(cm);
443			}
444
445	greg	2.15	/* Load and convert a matrix BSDF from the given XML file */
446			static CMATRIX *
447	greg	2.20	cm_loadBSDF(char *fname, COLOR cLamb)
448	greg	2.15	{
449	greg	2.21	CMATRIX *Tmat;
450			char *fpath;
451	greg	2.22	int recip;
452	greg	2.21	SDError ec;
453			SDData myBSDF;
454			SDSpectralDF *tdf;
455			COLOR bsdfLamb, specCol;
456	greg	2.15	/* find path to BSDF file */
457			fpath = getpath(fname, getrlibpath(), R_OK);
458			if (fpath == NULL) {
459			sprintf(errmsg, "cannot find BSDF file '%s'", fname);
460			error(USER, errmsg);
461			}
462	greg	2.16	SDclearBSDF(&myBSDF, fname); /* load XML and check type */
463	greg	2.15	ec = SDloadFile(&myBSDF, fpath);
464			if (ec)
465			error(USER, transSDError(ec));
466	greg	2.20	ccy2rgb(&myBSDF.tLamb.spec, myBSDF.tLamb.cieY/PI, bsdfLamb);
467	greg	2.23	recip = (myBSDF.tb == NULL);
468			tdf = recip ? myBSDF.tf : myBSDF.tb;
469	greg	2.21	if (tdf == NULL) { /* no non-Lambertian transmission? */
470	greg	2.20	if (cLamb != NULL)
471			copycolor(cLamb, bsdfLamb);
472			SDfreeBSDF(&myBSDF);
473			return(NULL);
474			}
475	greg	2.21	if (tdf->ncomp != 1 \|\| tdf->comp[0].func != &SDhandleMtx) {
476	greg	2.15	sprintf(errmsg, "unsupported BSDF '%s'", fpath);
477			error(USER, errmsg);
478			}
479			/* convert BTDF to matrix */
480	greg	2.21	ccy2rgb(&tdf->comp[0].cspec[0], 1., specCol);
481	greg	2.22	Tmat = recip ? cm_bsdf_recip(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist)
482			: cm_bsdf(bsdfLamb, specCol, (SDMat *)tdf->comp[0].dist);
483	greg	2.20	if (cLamb != NULL) /* Lambertian is included */
484			setcolor(cLamb, .0, .0, .0);
485	greg	2.15	/* free BSDF and return */
486			SDfreeBSDF(&myBSDF);
487			return(Tmat);
488			}
489
490	greg	2.24	/* Sum together a set of images and write result to fout */
491	greg	2.1	static int
492	greg	2.13	sum_images(const char fspec, const CMATRIX cv, FILE *fout)
493	greg	2.1	{
494			int myDT = DTfromHeader;
495	greg	2.24	COLOR *scanline = NULL;
496			CMATRIX *pmat = NULL;
497			int myXR=0, myYR=0;
498	greg	2.1	int i, y;
499
500			if (cv->ncols != 1)
501			error(INTERNAL, "expected vector in sum_images()");
502			for (i = 0; i < cv->nrows; i++) {
503			const COLORV *scv = cv_lval(cv,i);
504			char fname[1024];
505			FILE *fp;
506			int dt, xr, yr;
507			COLORV *psp;
508	greg	2.14	/* check for zero */
509	greg	2.19	if ((scv[RED] == 0) & (scv[GRN] == 0) & (scv[BLU] == 0) &&
510			(myDT != DTfromHeader) \| (i < cv->nrows-1))
511	greg	2.14	continue;
512	greg	2.1	/* open next picture */
513			sprintf(fname, fspec, i);
514			if ((fp = fopen(fname, "r")) == NULL) {
515			sprintf(errmsg, "cannot open picture '%s'", fname);
516			error(SYSTEM, errmsg);
517			}
518			SET_FILE_BINARY(fp);
519			dt = getDTfromHeader(fp);
520			if ((dt != DTrgbe) & (dt != DTxyze) \|\|
521			!fscnresolu(&xr, &yr, fp)) {
522			sprintf(errmsg, "file '%s' not a picture", fname);
523			error(USER, errmsg);
524			}
525			if (myDT == DTfromHeader) { /* on first one */
526			myDT = dt;
527			myXR = xr; myYR = yr;
528			scanline = (COLOR )malloc(sizeof(COLOR)myXR);
529			if (scanline == NULL)
530			error(SYSTEM, "out of memory in sum_images()");
531			pmat = cm_alloc(myYR, myXR);
532			memset(pmat->cmem, 0, sizeof(COLOR)myXRmyYR);
533	greg	2.2	/* finish header */
534	greg	2.13	fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, fout);
535			fputc('\n', fout);
536			fprtresolu(myXR, myYR, fout);
537			fflush(fout);
538	greg	2.1	} else if ((dt != myDT) \| (xr != myXR) \| (yr != myYR)) {
539			sprintf(errmsg, "picture '%s' format/size mismatch",
540			fname);
541			error(USER, errmsg);
542			}
543			psp = pmat->cmem;
544			for (y = 0; y < yr; y++) { /* read it in */
545			int x;
546			if (freadscan(scanline, xr, fp) < 0) {
547			sprintf(errmsg, "error reading picture '%s'",
548			fname);
549			error(SYSTEM, errmsg);
550			}
551			/* sum in scanline */
552			for (x = 0; x < xr; x++, psp += 3) {
553			multcolor(scanline[x], scv);
554			addcolor(psp, scanline[x]);
555			}
556			}
557			fclose(fp); /* done this picture */
558			}
559			free(scanline);
560			/* write scanlines */
561			for (y = 0; y < myYR; y++)
562	greg	2.13	if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0)
563	greg	2.1	return(0);
564			cm_free(pmat); /* all done */
565	greg	2.13	return(fflush(fout) == 0);
566	greg	2.1	}
567
568	greg	2.11	/* check to see if a string contains a %d or %o specification */
569			static int
570			hasNumberFormat(const char *s)
571	greg	2.1	{
572	greg	2.24	if (s == NULL)
573			return(0);
574
575	greg	2.18	while (*s) {
576			while (*s != '%')
577			if (!*s++)
578			return(0);
579			if (++s == '%') { / ignore "%%" */
580			++s;
581			continue;
582			}
583			while (isdigit(s)) / field length */
584			++s;
585			/* field we'll use? */
586			if ((s == 'd') \| (s == 'i') \| (*s == 'o') \|
587			(s == 'x') \| (s == 'X'))
588			return(1);
589			}
590			return(0); /* didn't find one */
591	greg	2.1	}
592
593			int
594			main(int argc, char *argv[])
595			{
596	greg	2.25	int skyfmt = DTascii;
597	greg	2.24	int nsteps = 1;
598			char *ofspec = NULL;
599			FILE *ofp = stdout;
600			CMATRIX cmtx; / component vector/matrix result */
601			char fnbuf[256];
602			int a, i;
603	greg	2.1
604			progname = argv[0];
605	greg	2.24	/* get options */
606	greg	2.25	for (a = 1; a < argc && argv[a][0] == '-'; a++)
607			switch (argv[a][1]) {
608	greg	2.24	case 'n':
609			nsteps = atoi(argv[++a]);
610			if (nsteps <= 0)
611			goto userr;
612			break;
613			case 'o':
614			ofspec = argv[++a];
615			break;
616	greg	2.25	case 'i':
617			switch (argv[a][2]) {
618			case 'f':
619			skyfmt = DTfloat;
620			break;
621			case 'd':
622			skyfmt = DTdouble;
623			break;
624			case 'a':
625			skyfmt = DTascii;
626			break;
627			default:
628			goto userr;
629			}
630			break;
631	greg	2.24	default:
632			goto userr;
633			}
634			if ((argc-a < 1) \| (argc-a > 4))
635			goto userr;
636	greg	2.1
637	greg	2.24	if (argc-a > 2) { /* VTDs expression */
638			CMATRIX smtx, Dmat, Tmat, imtx;
639	greg	2.20	COLOR tLamb;
640	greg	2.24	/* get sky vector/matrix */
641	greg	2.25	smtx = cm_load(argv[a+3], 0, nsteps, skyfmt);
642	greg	2.12	/* load BSDF */
643	greg	2.24	Tmat = cm_loadBSDF(argv[a+1], tLamb);
644	greg	2.5	/* load Daylight matrix */
645	greg	2.24	Dmat = cm_load(argv[a+2], Tmat==NULL ? 0 : Tmat->ncols,
646			smtx->nrows, DTfromHeader);
647	greg	2.1	/* multiply vector through */
648	greg	2.24	imtx = cm_multiply(Dmat, smtx);
649			cm_free(Dmat); cm_free(smtx);
650	greg	2.20	if (Tmat == NULL) { /* diffuse only */
651	greg	2.24	cmtx = cm_scale(imtx, tLamb);
652	greg	2.20	} else { /* else apply BTDF matrix */
653	greg	2.24	cmtx = cm_multiply(Tmat, imtx);
654	greg	2.20	cm_free(Tmat);
655			}
656	greg	2.24	cm_free(imtx);
657			} else { /* sky vector/matrix only */
658	greg	2.25	cmtx = cm_load(argv[a+1], 0, nsteps, skyfmt);
659	greg	2.24	}
660			/* prepare output stream */
661			if ((ofspec != NULL) & (nsteps == 1) && hasNumberFormat(ofspec)) {
662			sprintf(fnbuf, ofspec, 1);
663			ofspec = fnbuf;
664			}
665			if (ofspec != NULL && !hasNumberFormat(ofspec)) {
666			if ((ofp = fopen(ofspec, "w")) == NULL) {
667			fprintf(stderr, "%s: cannot open '%s' for output\n",
668			progname, ofspec);
669			return(1);
670			}
671			ofspec = NULL; /* only need to open once */
672	greg	2.12	}
673	greg	2.24	if (hasNumberFormat(argv[a])) { /* generating image(s) */
674			if (ofspec == NULL) {
675			SET_FILE_BINARY(ofp);
676			newheader("RADIANCE", ofp);
677			printargs(argc, argv, ofp);
678			fputnow(ofp);
679			}
680			if (nsteps > 1) /* multiple output frames? */
681			for (i = 0; i < nsteps; i++) {
682			CMATRIX *cvec = cm_column(cmtx, i);
683			if (ofspec != NULL) {
684			sprintf(fnbuf, ofspec, i+1);
685			if ((ofp = fopen(fnbuf, "wb")) == NULL) {
686			fprintf(stderr,
687			"%s: cannot open '%s' for output\n",
688			progname, fnbuf);
689			return(1);
690			}
691			newheader("RADIANCE", ofp);
692			printargs(argc, argv, ofp);
693			fputnow(ofp);
694			}
695			fprintf(ofp, "FRAME=%d\n", i+1);
696			if (!sum_images(argv[a], cvec, ofp))
697			return(1);
698			if (ofspec != NULL) {
699			if (fclose(ofp) == EOF) {
700			fprintf(stderr,
701			"%s: error writing to '%s'\n",
702			progname, fnbuf);
703			return(1);
704			}
705			ofp = stdout;
706			}
707			cm_free(cvec);
708			}
709			else if (!sum_images(argv[a], cmtx, ofp))
710	greg	2.1	return(1);
711	greg	2.24	} else { /* generating vector/matrix */
712			CMATRIX *Vmat = cm_load(argv[a], 0, cmtx->nrows, DTfromHeader);
713			CMATRIX *rmtx = cm_multiply(Vmat, cmtx);
714	greg	2.1	cm_free(Vmat);
715	greg	2.24	if (ofspec != NULL) /* multiple vector files? */
716			for (i = 0; i < nsteps; i++) {
717			CMATRIX *rvec = cm_column(rmtx, i);
718			sprintf(fnbuf, ofspec, i+1);
719			if ((ofp = fopen(fnbuf, "w")) == NULL) {
720			fprintf(stderr,
721			"%s: cannot open '%s' for output\n",
722			progname, fnbuf);
723			return(1);
724			}
725			cm_print(rvec, ofp);
726			if (fclose(ofp) == EOF) {
727			fprintf(stderr,
728			"%s: error writing to '%s'\n",
729			progname, fnbuf);
730			return(1);
731			}
732			ofp = stdout;
733			cm_free(rvec);
734			}
735			else
736			cm_print(rmtx, ofp);
737			cm_free(rmtx);
738			}
739			if (fflush(ofp) == EOF) { /* final clean-up */
740			fprintf(stderr, "%s: write error on output\n", progname);
741			return(1);
742	greg	2.1	}
743	greg	2.24	cm_free(cmtx);
744	greg	2.1	return(0);
745	greg	2.24	userr:
746	greg	2.25	fprintf(stderr, "Usage: %s [-n nsteps][-o ospec][-i{f\|d}] DCspec [skyf]\n",
747	greg	2.24	progname);
748	greg	2.25	fprintf(stderr, " or: %s [-n nsteps][-o ospec][-i{f\|d}] Vspec Tbsdf.xml Dmat.dat [skyf]\n",
749	greg	2.24	progname);
750			return(1);
751	greg	2.1	}