src/util/dctimestep.c

#ifndef lint
static const char RCSid[] = "$Id: dctimestep.c,v 2.8 2009/06/21 21:49:13 greg Exp $";
#endif
/*
 * Compute time-step result using Daylight Coefficient method.
 *
 *      G. Ward
 */

#include <ctype.h>
#include "standard.h"
#include "platform.h"
#include "color.h"
#include "resolu.h"
#include "bsdf.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))
                return(NULL);
        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)
{
        CMATRIX *cm;
        FILE    *fp = stdin;

        if (fname == NULL)
                fname = "<stdin>";
        else if ((fp = fopen(fname, "r")) == NULL) {
                sprintf(errmsg, "cannot open file '%s'", fname);
                error(SYSTEM, errmsg);
        }
        if (dtype != DTascii)
                SET_FILE_BINARY(fp);
        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)
                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 (getc(fp) != EOF) {
                                sprintf(errmsg,
                                "unexpected data at end of binary file %s",
                                                fname);
                                error(WARNING, errmsg);
                }
        }
        if (fp != stdin)
                fclose(fp);
        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 */
}

/* 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 struct BSDF_data *bsdf)
{
        CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
        int     nbadohm = 0;
        int     nneg = 0;
        int     r, c;
        
        for (c = 0; c < cm->ncols; c++) {
                float   dom = getBSDF_incohm(bsdf,c);
                FVECT   v;
                
                if (dom <= .0) {
                        nbadohm++;
                        continue;
                }
                if (!getBSDF_incvec(v,bsdf,c) || v[2] > FTINY)
                        error(USER, "illegal incoming BTDF direction");
                dom *= -v[2];

                for (r = 0; r < cm->nrows; r++) {
                        float   f = BSDF_value(bsdf,c,r);
                        COLORV  *mp = cm_lval(cm,r,c);

                        if (f <= .0) {
                                nneg += (f < -FTINY);
                                f = .0f;
                        }
                        mp[0] = mp[1] = mp[2] = f * dom;
                }
        }
        if (nneg || nbadohm) {
                sprintf(errmsg,
                    "BTDF has %d negatives and %d bad incoming solid angles",
                                nneg, nbadohm);
                error(WARNING, errmsg);
        }
        return(cm);
}

/* Sum together a set of images and write result to stdout */
static int
sum_images(const char *fspec, const CMATRIX *cv)
{
        int     myDT = DTfromHeader;
        CMATRIX *pmat;
        COLOR   *scanline;
        int     myXR, myYR;
        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;
                                                        /* 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, stdout);
                        fputc('\n', stdout);
                        fprtresolu(myXR, myYR, stdout);
                        fflush(stdout);
                } 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, stdout) < 0)
                        return(0);
        cm_free(pmat);                                  /* all done */
        return(fflush(stdout) == 0);
}

/* check to see if a string contains a %d specification */
int
hasDecimalSpec(const char *s)
{
        while (*s && *s != '%')
                s++;
        if (!*s)
                return(0);
        do
                ++s;
        while (isdigit(*s));

        return(*s == 'd');
}

int
main(int argc, char *argv[])
{
        CMATRIX                 *tvec, *Dmat, *Tmat, *ivec, *cvec;
        struct BSDF_data        *btdf;

        progname = argv[0];

        if ((argc < 4) | (argc > 5)) {
                fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
                                progname);
                return(1);
        }
                                                /* load Tregenza vector */
        tvec = cm_load(argv[4], 0, 1, DTascii); /* argv[4]==NULL iff argc==4 */
                                                /* load BTDF */
        btdf = load_BSDF(argv[2]);
        if (btdf == NULL)
                return(1);
                                                /* load Daylight matrix */
        Dmat = cm_load(argv[3], btdf->ninc, tvec->nrows, DTfromHeader);
                                                /* multiply vector through */
        ivec = 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, ivec);         /* cvec = component vector */
        cm_free(Tmat); cm_free(ivec);
        if (hasDecimalSpec(argv[1])) {          /* generating image */
                SET_FILE_BINARY(stdout);
                newheader("RADIANCE", stdout);
                printargs(argc, argv, stdout);
                fputnow(stdout);
                if (!sum_images(argv[1], cvec))
                        return(1);
        } else {                                /* generating vector */
                CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);
                CMATRIX *rvec = cm_multiply(Vmat, cvec);
                cm_free(Vmat);
                cm_print(rvec, stdout);
                cm_free(rvec);
        }
        cm_free(cvec);                          /* final clean-up */
        return(0);
}
Revision:	2.9
Committed:	Sun Jun 21 21:53:10 2009 UTC (16 years, 4 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.8:	+2 -3 lines
Log Message:	Minor tweaks
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.9	static const char RCSid[] = "$Id: dctimestep.c,v 2.8 2009/06/21 21:49:13 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			#include "color.h"
14			#include "resolu.h"
15			#include "bsdf.h"
16
17			char progname; / global argv[0] */
18
19	greg	2.2	/* Data types for file loading */
20	greg	2.1	enum {DTfromHeader, DTascii, DTfloat, DTdouble, DTrgbe, DTxyze};
21
22			/* A color coefficient matrix -- vectors have ncols==1 */
23			typedef struct {
24			int nrows, ncols;
25			COLORV cmem[3]; /* extends struct */
26			} CMATRIX;
27
28			#define COLSPEC (sizeof(COLORV)==sizeof(float) ? "%f %f %f" : "%lf %lf %lf")
29
30			#define cm_lval(cm,r,c) ((cm)->cmem + 3((r)(cm)->ncols + (c)))
31
32			#define cv_lval(cm,i) ((cm)->cmem + 3*(i))
33
34			/* Allocate a color coefficient matrix */
35			static CMATRIX *
36			cm_alloc(int nrows, int ncols)
37			{
38			CMATRIX *cm;
39
40			if ((nrows <= 0) \| (ncols <= 0))
41			return(NULL);
42			cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
43			3sizeof(COLORV)(nrows*ncols - 1));
44			if (cm == NULL)
45			error(SYSTEM, "out of memory in cm_alloc()");
46			cm->nrows = nrows;
47			cm->ncols = ncols;
48			return(cm);
49			}
50
51			#define cm_free(cm) free(cm)
52
53			/* Resize color coefficient matrix */
54			static CMATRIX *
55			cm_resize(CMATRIX *cm, int nrows)
56			{
57			if (nrows == cm->nrows)
58			return(cm);
59			if (nrows <= 0) {
60			cm_free(cm);
61			return(NULL);
62			}
63			cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) +
64			3sizeof(COLORV)(nrows*cm->ncols - 1));
65			if (cm == NULL)
66			error(SYSTEM, "out of memory in cm_resize()");
67			cm->nrows = nrows;
68			return(cm);
69			}
70
71			/* Load header to obtain data type */
72			static int
73			getDT(char s, void p)
74			{
75			char fmt[32];
76
77			if (formatval(fmt, s)) {
78			if (!strcmp(fmt, "ascii"))
79			((int )p) = DTascii;
80			else if (!strcmp(fmt, "float"))
81			((int )p) = DTfloat;
82			else if (!strcmp(fmt, "double"))
83			((int )p) = DTdouble;
84			else if (!strcmp(fmt, COLRFMT))
85			((int )p) = DTrgbe;
86			else if (!strcmp(fmt, CIEFMT))
87			((int )p) = DTxyze;
88			}
89			return(0);
90			}
91
92			static int
93			getDTfromHeader(FILE *fp)
94			{
95			int dt = DTfromHeader;
96
97			if (getheader(fp, getDT, &dt) < 0)
98			error(SYSTEM, "header read error");
99			if (dt == DTfromHeader)
100			error(USER, "missing data format in header");
101			return(dt);
102			}
103
104			/* Allocate and load a matrix from the given file (or stdin if NULL) */
105			static CMATRIX *
106			cm_load(const char *fname, int nrows, int ncols, int dtype)
107			{
108			CMATRIX *cm;
109			FILE *fp = stdin;
110
111			if (fname == NULL)
112			fname = "<stdin>";
113			else if ((fp = fopen(fname, "r")) == NULL) {
114			sprintf(errmsg, "cannot open file '%s'", fname);
115			error(SYSTEM, errmsg);
116			}
117			if (dtype != DTascii)
118			SET_FILE_BINARY(fp);
119			if (dtype == DTfromHeader)
120			dtype = getDTfromHeader(fp);
121			switch (dtype) {
122			case DTascii:
123			case DTfloat:
124			case DTdouble:
125			break;
126			default:
127			error(USER, "unexpected data type in cm_load()");
128			}
129			if (nrows <= 0) { /* don't know length? */
130			int guessrows = 147; /* usually big enough */
131			if ((dtype != DTascii) & (fp != stdin)) {
132			long startpos = ftell(fp);
133			if (fseek(fp, 0L, SEEK_END) == 0) {
134			long endpos = ftell(fp);
135			long elemsiz = 3*(dtype==DTfloat ?
136			sizeof(float) : sizeof(double));
137	greg	2.2
138			if ((endpos - startpos) % (ncols*elemsiz)) {
139			sprintf(errmsg,
140			"improper length for binary file '%s'",
141			fname);
142			error(USER, errmsg);
143			}
144			guessrows = (endpos - startpos)/(ncols*elemsiz);
145	greg	2.1	if (fseek(fp, startpos, SEEK_SET) < 0) {
146			sprintf(errmsg,
147			"fseek() error on file '%s'",
148			fname);
149			error(SYSTEM, errmsg);
150			}
151	greg	2.2	nrows = guessrows; /* we're confident */
152	greg	2.1	}
153			}
154			cm = cm_alloc(guessrows, ncols);
155			} else
156			cm = cm_alloc(nrows, ncols);
157			if (cm == NULL)
158			return(NULL);
159			if (dtype == DTascii) { /* read text file */
160			int maxrow = (nrows > 0 ? nrows : 32000);
161			int r, c;
162			for (r = 0; r < maxrow; r++) {
163	greg	2.6	if (r >= cm->nrows) /* need more space? */
164	greg	2.1	cm = cm_resize(cm, 2*cm->nrows);
165			for (c = 0; c < ncols; c++) {
166			COLORV *cv = cm_lval(cm,r,c);
167			if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3)
168	greg	2.6	if ((nrows <= 0) & (r > 0) & !c) {
169	greg	2.1	cm = cm_resize(cm, maxrow=r);
170			break;
171			} else
172			goto EOFerror;
173			}
174			}
175			while ((c = getc(fp)) != EOF)
176			if (!isspace(c)) {
177			sprintf(errmsg,
178			"unexpected data at end of ascii file %s",
179			fname);
180			error(WARNING, errmsg);
181			break;
182			}
183			} else { /* read binary file */
184			if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) :
185			sizeof(double))) {
186			int nread = 0;
187			do { /* read all we can */
188			nread += fread(cm->cmem + 3*nread,
189			3*sizeof(COLORV),
190			cm->nrows*cm->ncols - nread,
191			fp);
192			if (nrows <= 0) { /* unknown length */
193			if (nread == cm->nrows*cm->ncols)
194			/* need more space? */
195			cm = cm_resize(cm, 2*cm->nrows);
196	greg	2.6	else if (nread && !(nread % cm->ncols))
197	greg	2.1	/* seem to be done */
198			cm = cm_resize(cm, nread/cm->ncols);
199			else /* ended mid-row */
200			goto EOFerror;
201			} else if (nread < cm->nrows*cm->ncols)
202			goto EOFerror;
203			} while (nread < cm->nrows*cm->ncols);
204
205			} else if (dtype == DTdouble) {
206			double dc[3]; /* load from double */
207			COLORV *cvp = cm->cmem;
208			int n = nrows*ncols;
209
210			if (n <= 0)
211			goto not_handled;
212			while (n--) {
213			if (fread(dc, sizeof(double), 3, fp) != 3)
214			goto EOFerror;
215			copycolor(cvp, dc);
216			cvp += 3;
217			}
218			} else /* dtype == DTfloat */ {
219			float fc[3]; /* load from float */
220			COLORV *cvp = cm->cmem;
221			int n = nrows*ncols;
222
223			if (n <= 0)
224			goto not_handled;
225			while (n--) {
226			if (fread(fc, sizeof(float), 3, fp) != 3)
227			goto EOFerror;
228			copycolor(cvp, fc);
229			cvp += 3;
230			}
231			}
232			if (getc(fp) != EOF) {
233			sprintf(errmsg,
234			"unexpected data at end of binary file %s",
235			fname);
236			error(WARNING, errmsg);
237			}
238			}
239			if (fp != stdin)
240			fclose(fp);
241			return(cm);
242			EOFerror:
243			sprintf(errmsg, "unexpected EOF reading %s",
244			fname);
245			error(USER, errmsg);
246			not_handled:
247			error(INTERNAL, "unhandled data size or length in cm_load()");
248			return(NULL); /* gratis return */
249			}
250
251			/* Multiply two matrices (or a matrix and a vector) and allocate the result*/
252			static CMATRIX *
253			cm_multiply(const CMATRIX cm1, const CMATRIX cm2)
254			{
255			CMATRIX *cmr;
256			int dr, dc, i;
257
258			if ((cm1->ncols <= 0) \| (cm1->ncols != cm2->nrows))
259			error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
260			cmr = cm_alloc(cm1->nrows, cm2->ncols);
261			if (cmr == NULL)
262			return(NULL);
263			for (dr = 0; dr < cmr->nrows; dr++)
264			for (dc = 0; dc < cmr->ncols; dc++) {
265			COLORV *dp = cm_lval(cmr,dr,dc);
266			dp[0] = dp[1] = dp[2] = 0;
267			for (i = 0; i < cm1->ncols; i++) {
268			const COLORV *cp1 = cm_lval(cm1,dr,i);
269			const COLORV *cp2 = cm_lval(cm2,i,dc);
270			dp[0] += cp1[0] * cp2[0];
271			dp[1] += cp1[1] * cp2[1];
272			dp[2] += cp1[2] * cp2[2];
273			}
274			}
275			return(cmr);
276			}
277
278			/* print out matrix as ASCII text -- no header */
279			static void
280			cm_print(const CMATRIX cm, FILE fp)
281			{
282			int r, c;
283			const COLORV *mp = cm->cmem;
284
285			for (r = 0; r < cm->nrows; r++) {
286			for (c = 0; c < cm->ncols; c++, mp += 3)
287			fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]);
288			fputc('\n', fp);
289			}
290			}
291
292			/* convert a BSDF to our matrix representation */
293			static CMATRIX *
294			cm_bsdf(const struct BSDF_data *bsdf)
295			{
296			CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
297	greg	2.3	int nbadohm = 0;
298			int nneg = 0;
299	greg	2.1	int r, c;
300
301	greg	2.8	for (c = 0; c < cm->ncols; c++) {
302	greg	2.7	float dom = getBSDF_incohm(bsdf,c);
303			FVECT v;
304
305			if (dom <= .0) {
306			nbadohm++;
307			continue;
308			}
309			if (!getBSDF_incvec(v,bsdf,c) \|\| v[2] > FTINY)
310			error(USER, "illegal incoming BTDF direction");
311			dom *= -v[2];
312
313			for (r = 0; r < cm->nrows; r++) {
314	greg	2.3	float f = BSDF_value(bsdf,c,r);
315	greg	2.9	COLORV *mp = cm_lval(cm,r,c);
316	greg	2.7
317	greg	2.3	if (f <= .0) {
318			nneg += (f < -FTINY);
319	greg	2.8	f = .0f;
320	greg	2.3	}
321	greg	2.7	mp[0] = mp[1] = mp[2] = f * dom;
322	greg	2.3	}
323	greg	2.7	}
324	greg	2.4	if (nneg \|\| nbadohm) {
325	greg	2.3	sprintf(errmsg,
326	greg	2.4	"BTDF has %d negatives and %d bad incoming solid angles",
327			nneg, nbadohm);
328	greg	2.3	error(WARNING, errmsg);
329			}
330	greg	2.1	return(cm);
331			}
332
333			/* Sum together a set of images and write result to stdout */
334			static int
335			sum_images(const char fspec, const CMATRIX cv)
336			{
337			int myDT = DTfromHeader;
338			CMATRIX *pmat;
339			COLOR *scanline;
340			int myXR, myYR;
341			int i, y;
342
343			if (cv->ncols != 1)
344			error(INTERNAL, "expected vector in sum_images()");
345			for (i = 0; i < cv->nrows; i++) {
346			const COLORV *scv = cv_lval(cv,i);
347			char fname[1024];
348			FILE *fp;
349			int dt, xr, yr;
350			COLORV *psp;
351			/* open next picture */
352			sprintf(fname, fspec, i);
353			if ((fp = fopen(fname, "r")) == NULL) {
354			sprintf(errmsg, "cannot open picture '%s'", fname);
355			error(SYSTEM, errmsg);
356			}
357			SET_FILE_BINARY(fp);
358			dt = getDTfromHeader(fp);
359			if ((dt != DTrgbe) & (dt != DTxyze) \|\|
360			!fscnresolu(&xr, &yr, fp)) {
361			sprintf(errmsg, "file '%s' not a picture", fname);
362			error(USER, errmsg);
363			}
364			if (myDT == DTfromHeader) { /* on first one */
365			myDT = dt;
366			myXR = xr; myYR = yr;
367			scanline = (COLOR )malloc(sizeof(COLOR)myXR);
368			if (scanline == NULL)
369			error(SYSTEM, "out of memory in sum_images()");
370			pmat = cm_alloc(myYR, myXR);
371			memset(pmat->cmem, 0, sizeof(COLOR)myXRmyYR);
372	greg	2.2	/* finish header */
373			fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, stdout);
374			fputc('\n', stdout);
375			fprtresolu(myXR, myYR, stdout);
376			fflush(stdout);
377	greg	2.1	} else if ((dt != myDT) \| (xr != myXR) \| (yr != myYR)) {
378			sprintf(errmsg, "picture '%s' format/size mismatch",
379			fname);
380			error(USER, errmsg);
381			}
382			psp = pmat->cmem;
383			for (y = 0; y < yr; y++) { /* read it in */
384			int x;
385			if (freadscan(scanline, xr, fp) < 0) {
386			sprintf(errmsg, "error reading picture '%s'",
387			fname);
388			error(SYSTEM, errmsg);
389			}
390			/* sum in scanline */
391			for (x = 0; x < xr; x++, psp += 3) {
392			multcolor(scanline[x], scv);
393			addcolor(psp, scanline[x]);
394			}
395			}
396			fclose(fp); /* done this picture */
397			}
398			free(scanline);
399			/* write scanlines */
400			for (y = 0; y < myYR; y++)
401			if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, stdout) < 0)
402			return(0);
403			cm_free(pmat); /* all done */
404			return(fflush(stdout) == 0);
405			}
406
407			/* check to see if a string contains a %d specification */
408			int
409			hasDecimalSpec(const char *s)
410			{
411			while (s && s != '%')
412			s++;
413			if (!*s)
414			return(0);
415			do
416			++s;
417			while (isdigit(*s));
418
419			return(*s == 'd');
420			}
421
422			int
423			main(int argc, char *argv[])
424			{
425	greg	2.2	CMATRIX tvec, Dmat, Tmat, ivec, *cvec;
426	greg	2.1	struct BSDF_data *btdf;
427
428			progname = argv[0];
429
430			if ((argc < 4) \| (argc > 5)) {
431			fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
432			progname);
433			return(1);
434			}
435	greg	2.5	/* load Tregenza vector */
436	greg	2.1	tvec = cm_load(argv[4], 0, 1, DTascii); /* argv[4]==NULL iff argc==4 */
437	greg	2.5	/* load BTDF */
438	greg	2.1	btdf = load_BSDF(argv[2]);
439			if (btdf == NULL)
440			return(1);
441	greg	2.5	/* load Daylight matrix */
442			Dmat = cm_load(argv[3], btdf->ninc, tvec->nrows, DTfromHeader);
443	greg	2.1	/* multiply vector through */
444	greg	2.2	ivec = cm_multiply(Dmat, tvec);
445	greg	2.1	cm_free(Dmat); cm_free(tvec);
446			Tmat = cm_bsdf(btdf); /* convert BTDF to matrix */
447			free_BSDF(btdf);
448	greg	2.2	cvec = cm_multiply(Tmat, ivec); /* cvec = component vector */
449			cm_free(Tmat); cm_free(ivec);
450	greg	2.1	if (hasDecimalSpec(argv[1])) { /* generating image */
451			SET_FILE_BINARY(stdout);
452			newheader("RADIANCE", stdout);
453			printargs(argc, argv, stdout);
454			fputnow(stdout);
455			if (!sum_images(argv[1], cvec))
456			return(1);
457			} else { /* generating vector */
458			CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);
459			CMATRIX *rvec = cm_multiply(Vmat, cvec);
460			cm_free(Vmat);
461			cm_print(rvec, stdout);
462			cm_free(rvec);
463			}
464			cm_free(cvec); /* final clean-up */
465			return(0);
466			}