ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxcomb.c
Revision: 2.6
Committed: Mon Dec 11 15:13:39 2023 UTC (5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +4 -4 lines
Log Message:
perf(rmtxcomb,pcomb): Avoid calling varset() if unnecessary

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.6 static const char RCSid[] = "$Id: rmtxcomb.c,v 2.5 2023/12/08 00:12:31 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * General component matrix combiner, operating on a row at a time.
6     */
7    
8     #include <errno.h>
9     #include <math.h>
10     #include "platform.h"
11     #include "rtio.h"
12     #include "resolu.h"
13     #include "rmatrix.h"
14     #include "calcomp.h"
15     #include "paths.h"
16    
17     #ifndef M_PI
18     #define M_PI 3.14159265358979323846
19     #endif
20    
21     #define MAXCOMP MAXCSAMP /* #components we support */
22    
23     /* Unary matrix operation(s) */
24     typedef struct {
25     double cmat[MAXCOMP*MAXCOMP]; /* component transformation */
26     double sca[MAXCOMP]; /* scalar coefficients */
27     const char *csym; /* symbolic coefficients */
28     short clen; /* number of coefficients */
29     short nsf; /* number of scalars */
30     } RUNARYOP;
31    
32     /* Input matrix */
33     typedef struct {
34     const char *inspec; /* input specification */
35     RUNARYOP preop; /* transform operation */
36     RMATRIX imx; /* input matrix header info */
37     RMATRIX *rmp; /* active single-row matrix */
38     FILE *infp; /* open input stream */
39     } ROPMAT;
40    
41     ROPMAT *mop = NULL; /* allocated input array */
42     int nall = 0; /* number allocated */
43     int nmats = 0; /* number of actual inputs */
44    
45     RMATRIX *mcat = NULL; /* final concatenation */
46     int mcat_last = 0; /* goes after trailing ops? */
47    
48     int in_nrows; /* input row count */
49     #define in_ncols (mop[0].rmp->ncols) /* input column count */
50     #define in_ncomp (mop[0].rmp->ncomp) /* input #components */
51    
52     extern int nowarn; /* turn off warnings? */
53    
54     int cur_row; /* current input/output row */
55     int cur_col; /* current input/output column */
56     int cur_chan; /* if we're looping channels */
57    
58     static int checksymbolic(ROPMAT *rop);
59    
60     static int
61     split_input(ROPMAT *rop)
62     {
63     if (rop->rmp == &rop->imx && !(rop->rmp = rmx_copy(&rop->imx))) {
64     fputs("Out of memory in split_input()\n", stderr);
65     return(0);
66     }
67     rmx_reset(rop->rmp);
68     return(1);
69     }
70    
71     /* Check/set transform based on a reference input file */
72     static int
73     checkreffile(ROPMAT *rop)
74     {
75     static const char *curRF = NULL;
76     static RMATRIX refm;
77     const int nc = rop->imx.ncomp;
78     int i;
79    
80     if (!curRF || strcmp(rop->preop.csym, curRF)) {
81     FILE *fp = fopen(rop->preop.csym, "rb");
82     if (!rmx_load_header(&refm, fp)) {
83     fprintf(stderr, "%s: cannot read info header\n",
84     rop->preop.csym);
85     curRF = NULL;
86     if (fp) fclose(fp);
87     return(0);
88     }
89     fclose(fp);
90     curRF = rop->preop.csym;
91     }
92 greg 2.5 if (refm.ncomp == 3) {
93 greg 2.1 rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
94     return(checksymbolic(rop));
95     }
96     if (refm.ncomp == 2) {
97     fprintf(stderr, "%s: cannot convert to 2 components\n",
98     curRF);
99     return(0);
100     }
101     if (refm.ncomp == 1) {
102     rop->preop.csym = "Y"; /* XXX big assumption */
103     return(checksymbolic(rop));
104     }
105     if (refm.ncomp == nc &&
106     !memcmp(refm.wlpart, rop->imx.wlpart, sizeof(refm.wlpart)))
107     return(1); /* nothing to do */
108    
109     if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
110     fprintf(stderr, "%s: cannot resample from %d to %d components\n",
111     curRF, nc, refm.ncomp);
112     return(0);
113     }
114     if (!split_input(rop)) /* get our own struct */
115     return(0);
116     rop->preop.clen = refm.ncomp * nc; /* compute spec to ref */
117    
118     for (i = 0; i < nc; i++) {
119     SCOLOR scstim, scresp;
120     int j;
121     memset(scstim, 0, sizeof(COLORV)*nc);
122     scstim[i] = 1.f;
123     convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
124     scstim, nc, rop->imx.wlpart[0], rop->imx.wlpart[3]);
125     for (j = refm.ncomp; j-- > 0; )
126     rop->preop.cmat[j*nc + i] = scresp[j];
127     }
128     /* remember new spectral params */
129     memcpy(rop->rmp->wlpart, refm.wlpart, sizeof(rop->rmp->wlpart));
130     rop->rmp->ncomp = refm.ncomp;
131     return(1);
132     }
133    
134     /* Compute conversion row from spectrum to one channel of RGB */
135     static void
136     rgbrow(ROPMAT *rop, int r, int p)
137     {
138     const int nc = rop->imx.ncomp;
139     const float * wlp = rop->imx.wlpart;
140     int i;
141    
142     for (i = nc; i--; ) {
143     int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
144     int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
145     COLOR crgb;
146     spec_rgb(crgb, nmStart, nmEnd);
147     rop->preop.cmat[r*nc+i] = crgb[p];
148     }
149     }
150    
151     /* Compute conversion row from spectrum to one channel of XYZ */
152     static void
153     xyzrow(ROPMAT *rop, int r, int p)
154     {
155     const int nc = rop->imx.ncomp;
156     const float * wlp = rop->imx.wlpart;
157     int i;
158    
159     for (i = nc; i--; ) {
160     int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
161     int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
162     COLOR cxyz;
163     spec_cie(cxyz, nmStart, nmEnd);
164     rop->preop.cmat[r*nc+i] = cxyz[p];
165     }
166     }
167    
168     /* Use the spectral sensitivity function to compute matrix coefficients */
169     static void
170     sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
171     {
172     const int nc = rop->imx.ncomp;
173     int i;
174    
175     for (i = nc; i--; ) {
176     SCOLOR sclr;
177     memset(sclr, 0, sizeof(COLORV)*nc);
178     sclr[i] = 1.f;
179     rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->imx.wlpart);
180     }
181     }
182    
183     /* Check/set symbolic transform */
184     static int
185     checksymbolic(ROPMAT *rop)
186     {
187     const int nc = rop->imx.ncomp;
188     const int dt = rop->imx.dtype;
189     int i, j;
190     /* check suffix => reference file */
191     if (strchr(rop->preop.csym, '.') > rop->preop.csym)
192     return(checkreffile(rop));
193    
194     if (nc < 3) {
195     fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
196     rop->inspec, rop->preop.csym);
197     return(0);
198     }
199     rop->preop.clen = strlen(rop->preop.csym) * nc;
200     if (rop->preop.clen > MAXCOMP*MAXCOMP) {
201     fprintf(stderr, "%s: -c '%s' results in too many components\n",
202     rop->inspec, rop->preop.csym);
203     return(0);
204     }
205     for (j = 0; rop->preop.csym[j]; j++) {
206     int comp = 0;
207     switch (rop->preop.csym[j]) {
208     case 'B':
209     ++comp;
210     /* fall through */
211     case 'G':
212     ++comp;
213     /* fall through */
214     case 'R':
215     if (dt == DTxyze) {
216     for (i = 3; i--; )
217     rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
218     xyz2rgbmat[comp][i];
219     } else if (nc == 3)
220     rop->preop.cmat[j*nc+comp] = 1.;
221     else
222     rgbrow(rop, j, comp);
223     break;
224     case 'Z':
225     ++comp;
226     /* fall through */
227     case 'Y':
228     ++comp;
229     /* fall through */
230     case 'X':
231     if (dt == DTxyze) {
232     rop->preop.cmat[j*nc+comp] = 1.;
233     } else if (nc == 3) {
234     for (i = 3; i--; )
235     rop->preop.cmat[j*nc+i] =
236     rgb2xyzmat[comp][i];
237     } else if (comp == CIEY)
238     sensrow(rop, j, scolor2photopic);
239     else
240     xyzrow(rop, j, comp);
241    
242     for (i = nc*(dt != DTxyze); i--; )
243     rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
244     break;
245     case 'S': /* scotopic (il)luminance */
246     sensrow(rop, j, scolor2scotopic);
247     for (i = nc; i--; )
248     rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
249     break;
250     case 'M': /* melanopic (il)luminance */
251     sensrow(rop, j, scolor2melanopic);
252     for (i = nc; i--; )
253     rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
254     break;
255     case 'A': /* average component */
256     for (i = nc; i--; )
257     rop->preop.cmat[j*nc+i] = 1./(double)nc;
258     break;
259     default:
260     fprintf(stderr, "%s: -c '%c' unsupported\n",
261     rop->inspec, rop->preop.csym[j]);
262     return(0);
263     }
264     }
265     if (!split_input(rop)) /* get our own struct */
266     return(0);
267     memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
268     rop->rmp->ncomp = rop->preop.clen / nc;
269     /* decide on output type */
270     if (!strcmp(rop->preop.csym, "XYZ")) {
271     if (dt <= DTspec)
272     rop->rmp->dtype = DTxyze;
273     } else if (!strcmp(rop->preop.csym, "RGB")) {
274     if (dt <= DTspec)
275     rop->rmp->dtype = DTrgbe;
276 greg 2.2 } else if (rop->rmp->dtype == DTspec)
277 greg 2.1 rop->rmp->dtype = DTfloat;
278     return(1);
279     }
280    
281     static int
282     get_component_xfm(ROPMAT *rop)
283     {
284     int i, j;
285    
286     if (rop->rmp != &rop->imx) { /* reset destination matrix */
287     rmx_free(rop->rmp);
288     rop->rmp = &rop->imx;
289     }
290     if (rop->preop.csym && /* symbolic transform? */
291     !checksymbolic(rop))
292     return(0);
293     /* undo exposure? */
294     if (fabs(1. - bright(rop->rmp->cexp)) > .025) {
295     if (rop->rmp->ncomp == 1)
296     rop->rmp->cexp[RED] = rop->rmp->cexp[GRN] =
297     rop->rmp->cexp[BLU] = bright(rop->rmp->cexp);
298     if (rop->preop.nsf <= 0) {
299     rop->preop.nsf = i = rop->rmp->ncomp;
300     while (i--)
301     rop->preop.sca[i] = 1.;
302     }
303 greg 2.2 if (rop->preop.nsf == 1) {
304     if (rop->rmp->ncomp == 3) {
305     rop->preop.sca[2] = rop->preop.sca[1] =
306     rop->preop.sca[0];
307     rop->preop.nsf = 3;
308     } else
309     rop->preop.sca[0] /= bright(rop->rmp->cexp);
310     }
311     if (rop->preop.nsf == 3) {
312 greg 2.1 opcolor(rop->preop.sca, /=, rop->rmp->cexp);
313 greg 2.2 } else if (rop->preop.nsf > 3) { /* punt */
314 greg 2.1 double mult = 1./bright(rop->rmp->cexp);
315     for (i = rop->preop.nsf; i--; )
316     rop->preop.sca[i] *= mult;
317     }
318     setcolor(rop->rmp->cexp, 1., 1., 1.);
319     }
320     if (rop->preop.clen > 0) { /* use component transform? */
321     if (rop->preop.clen % rop->imx.ncomp) {
322     fprintf(stderr, "%s: -c must have N x %d coefficients\n",
323     rop->inspec, rop->imx.ncomp);
324     return(0);
325     }
326     if (rop->preop.nsf > 0) { /* scale transform, instead */
327     if (rop->preop.nsf == 1) {
328     for (i = rop->preop.clen; i--; )
329     rop->preop.cmat[i] *= rop->preop.sca[0];
330     } else if (rop->preop.nsf*rop->imx.ncomp != rop->preop.clen) {
331     fprintf(stderr, "%s: -s must have one or %d factors\n",
332     rop->inspec,
333     rop->preop.clen/rop->imx.ncomp);
334     return(0);
335     } else {
336     for (i = rop->preop.nsf; i--; )
337     for (j = rop->imx.ncomp; j--; )
338     rop->preop.cmat[i*rop->imx.ncomp+j]
339     *= rop->preop.sca[i];
340     }
341     }
342     rop->preop.nsf = 0; /* now folded in */
343     if (!split_input(rop)) /* get our own struct */
344     return(0);
345     rop->rmp->ncomp = rop->preop.clen / rop->imx.ncomp;
346 greg 2.4 if ((rop->rmp->ncomp > 3) & (rop->rmp->dtype <= DTspec)) {
347 greg 2.1 rop->rmp->dtype = DTfloat; /* probably not actual spectrum */
348 greg 2.4 memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
349     }
350 greg 2.1 } else if (rop->preop.nsf > 0) { /* else use scalar(s)? */
351     if (rop->preop.nsf == 1) {
352     for (i = rop->rmp->ncomp; --i; )
353     rop->preop.sca[i] = rop->preop.sca[0];
354     rop->preop.nsf = rop->rmp->ncomp;
355     } else if (rop->preop.nsf != rop->rmp->ncomp) {
356     fprintf(stderr, "%s: -s must have one or %d factors\n",
357     rop->inspec, rop->rmp->ncomp);
358     return(0);
359     }
360     }
361     return(1);
362     }
363    
364     static int
365     apply_op(RMATRIX *dst, const RMATRIX *src, const RUNARYOP *ro)
366     {
367     if (ro->clen > 0) {
368     RMATRIX *res = rmx_transform(src, dst->ncomp, ro->cmat);
369     if (!res) {
370     fputs("Error in call to rmx_transform()\n", stderr);
371     return(0);
372     }
373 greg 2.4 if (!rmx_transfer_data(dst, res, 0))
374     return(0);
375 greg 2.1 rmx_free(res);
376     } else if (dst != src)
377     memcpy(dst->mtx, src->mtx,
378     sizeof(double)*dst->ncomp*dst->ncols*dst->nrows);
379     if (ro->nsf == dst->ncomp)
380     rmx_scale(dst, ro->sca);
381     return(1);
382     }
383    
384     static int
385     open_input(ROPMAT *rop)
386     {
387     int outtype;
388    
389     if (!rop || !rop->inspec || !rop->inspec[0])
390     return(0);
391     if (rop->inspec == stdin_name)
392     rop->infp = stdin;
393     else if (rop->inspec[0] == '!')
394     rop->infp = popen(rop->inspec+1, "r");
395     else
396     rop->infp = fopen(rop->inspec, "rb");
397    
398     if (!rmx_load_header(&rop->imx, rop->infp)) {
399     fprintf(stderr, "Bad header from: %s\n", rop->inspec);
400     return(0);
401     }
402     return(get_component_xfm(rop));
403     }
404    
405     /* Return nominal wavelength associated with input component (return nm) */
406     static double
407     l_wavelength(char *nam)
408     {
409     double comp = argument(1);
410    
411     if ((comp < -.5) | (comp >= in_ncomp+.5)) {
412     errno = EDOM;
413     return(.0);
414     }
415     if (comp < .5) /* asking for #components? */
416     return(in_ncomp);
417    
418     if (in_ncomp == 3) { /* special case for RGB */
419     const int w0 = (int)(comp - .5);
420     return(mop[0].rmp->wlpart[w0] +
421     (comp-.5)*(mop[0].rmp->wlpart[w0+1] -
422     mop[0].rmp->wlpart[w0]));
423     }
424     return(mop[0].rmp->wlpart[0] + /* general case, even div. */
425     (comp-.5)/(double)in_ncomp *
426     (mop[0].rmp->wlpart[3] - mop[0].rmp->wlpart[0]));
427     }
428    
429     /* Return ith input with optional channel selector */
430     static double
431     l_chanin(char *nam)
432     {
433     double inp = argument(1);
434     int mi, chan;
435    
436     if ((mi = (int)(inp-.5)) < 0 || mi >= nmats) {
437     errno = EDOM;
438     return(.0);
439     }
440 greg 2.2 if (inp < .5) /* asking for #inputs? */
441 greg 2.1 return(nmats);
442    
443     if (nargum() >= 2) {
444     double cval = argument(2);
445     if (cval < .5 || (chan = (int)(cval-.5)) >= in_ncomp) {
446     errno = EDOM;
447     return(.0);
448     }
449     } else
450     chan = cur_chan;
451    
452     return(mop[mi].rmp->mtx[cur_col*in_ncomp + chan]);
453     }
454    
455     static int
456     initialize(RMATRIX *imp)
457     {
458     int i;
459     /* XXX struct is zeroed coming in */
460     setcolor(imp->cexp, 1.f, 1.f, 1.f);
461     for (i = 0; i < nmats; i++) { /* open each input */
462     int restype;
463     if (!open_input(&mop[i]))
464     return(0);
465     restype = mop[i].rmp->dtype;
466     if (!imp->dtype || (restype = rmx_newtype(restype, imp->dtype)) > 0)
467     imp->dtype = restype;
468     else
469     fprintf(stderr, "%s: warning - data type mismatch\n",
470     mop[i].inspec);
471     if (!i) {
472     imp->nrows = in_nrows = mop[0].rmp->nrows;
473     imp->ncols = mop[0].rmp->ncols;
474     imp->ncomp = mop[0].rmp->ncomp;
475     memcpy(imp->wlpart, mop[0].rmp->wlpart, sizeof(imp->wlpart));
476     } else if ((mop[i].rmp->nrows != imp->nrows) |
477     (mop[i].rmp->ncols != imp->ncols) |
478     (mop[i].rmp->ncomp != imp->ncomp)) {
479     fprintf(stderr, "%s: mismatch in size or #components\n",
480     mop[i].inspec);
481     return(0);
482     } /* XXX should check wlpart? */
483     } /* set up .cal environment */
484     esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
485     esupport &= ~(E_OUTCHAN|E_INCHAN);
486     varset("PI", ':', M_PI);
487     varset("nfiles", ':', nmats);
488     varset("nrows", ':', in_nrows);
489     varset("ncols", ':', in_ncols);
490     varset("ncomp", ':', in_ncomp);
491     varset("R", ':', 1.);
492     varset("G", ':', 2.);
493     varset("B", ':', 3.);
494     funset("wl", 1, ':', l_wavelength);
495     funset("ci", 1, '=', l_chanin);
496     scompile("ri(i)=ci(i,R);gi(i)=ci(i,G);bi(i)=ci(i,B)", NULL, 0);
497     return(1);
498     }
499    
500     static void
501     output_headinfo(FILE *fp)
502     {
503     int i;
504    
505     for (i = 0; i < nmats; i++) {
506     const char *cp = mop[i].imx.info;
507     fputs(mop[i].inspec, fp);
508     fputs(":\n", fp);
509     if (!cp) continue;
510     while (*cp) {
511     if (*cp == '\n') {
512     cp++; /* avoid inadvertant terminus */
513     continue;
514     }
515     fputc('\t', fp); /* indent this input's info */
516     do
517     putc(*cp, fp);
518     while (*cp++ != '\n');
519     }
520     }
521     }
522    
523     static int
524     combine_input(ROPMAT *res, FILE *fout)
525     {
526 greg 2.2 int set_r, set_c;
527 greg 2.1 RMATRIX *tmp = NULL;
528     int co_set;
529     int i;
530     /* allocate input row buffers */
531     for (i = 0; i < nmats; i++) {
532     mop[i].imx.nrows = 1; /* we'll be doing a row at a time */
533     if (!rmx_prepare(&mop[i].imx))
534     goto memerror;
535     if (mop[i].rmp != &mop[i].imx) {
536     mop[i].rmp->nrows = 1;
537     if (!rmx_prepare(mop[i].rmp))
538     goto memerror;
539     }
540     }
541     /* prep output row buffers */
542     if (mcat || res->preop.clen > 0) {
543     if (!split_input(res)) /* need separate buffer */
544     return(0);
545     if (res->preop.clen > 0)
546     res->rmp->ncomp = res->preop.clen / res->imx.ncomp;
547     res->rmp->nrows = 1;
548     if (!mcat | !mcat_last && !rmx_prepare(res->rmp))
549     goto memerror;
550     }
551     if (mcat && mcat_last &&
552 greg 2.3 !(tmp = rmx_alloc(1, res->imx.ncols, res->rmp->ncomp)))
553 greg 2.1 goto memerror;
554     res->imx.nrows = 1;
555     if (!rmx_prepare(&res->imx))
556     goto memerror;
557     /* figure out what the user set */
558     co_set = fundefined("co");
559     if (!co_set)
560     co_set = -vardefined("co");
561     if (!co_set & (in_ncomp == 3) && vardefined("ro") &&
562     vardefined("go") && vardefined("bo")) {
563     scompile("co(p)=select(p,ro,go,bo)", NULL, 0);
564     co_set = 1;
565     }
566 greg 2.6 if (co_set) { /* set if user wants, didn't set */
567     set_r = varlookup("r") != NULL && !vardefined("r");
568     set_c = varlookup("c") != NULL && !vardefined("c");
569 greg 2.2 } else /* save a little time */
570     set_r = set_c = 0;
571 greg 2.1 /* read/process row-by-row */
572     for (cur_row = 0; cur_row < in_nrows; cur_row++) {
573     RMATRIX *mres = NULL;
574     for (i = 0; i < nmats; i++) {
575 greg 2.2 if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
576     fprintf(stderr, "%s: read error at row %d\n",
577     mop[i].inspec, cur_row);
578 greg 2.1 return(0);
579 greg 2.2 }
580 greg 2.1 if (!apply_op(mop[i].rmp, &mop[i].imx, &mop[i].preop))
581     return(0);
582     }
583 greg 2.2 if (set_r) varset("r", '=', cur_row);
584 greg 2.1 for (cur_col = 0; cur_col < in_ncols; cur_col++) {
585 greg 2.2 if (set_c) varset("c", '=', cur_col);
586 greg 2.1 for (cur_chan = 0; cur_chan < in_ncomp; cur_chan++) {
587     const int ndx = cur_col*in_ncomp + cur_chan;
588     eclock++;
589     if (!co_set) { /* just summing elements? */
590     res->imx.mtx[ndx] = 0;
591     for (i = nmats; i--; )
592     res->imx.mtx[ndx] += mop[i].rmp->mtx[ndx];
593     } else if (co_set > 0) {
594 greg 2.2 double dchan = cur_chan+1;
595 greg 2.1 res->imx.mtx[ndx] = funvalue("co", 1, &dchan);
596     } else
597     res->imx.mtx[ndx] = varvalue("co");
598     }
599     } /* final conversions */
600     if (!mcat) {
601     if (!apply_op(res->rmp, &res->imx, &res->preop))
602     return(0);
603     } else if (mcat_last) {
604     if (!apply_op(tmp, &res->imx, &res->preop))
605     return(0);
606     mres = rmx_multiply(tmp, mcat);
607     if (!mres)
608     goto multerror;
609 greg 2.4 if (!rmx_transfer_data(res->rmp, mres, 0))
610     return(0);
611 greg 2.1 } else /* mcat && !mcat_last */ {
612     mres = rmx_multiply(&res->imx, mcat);
613     if (!mres)
614     goto multerror;
615     if (!apply_op(res->rmp, mres, &res->preop))
616     return(0);
617     }
618     rmx_free(mres); mres = NULL;
619     if (!rmx_write_data(res->rmp->mtx, res->rmp->ncomp,
620     res->rmp->ncols, res->rmp->dtype, fout))
621     return(0);
622     }
623 greg 2.2 #if 0 /* we're about to exit, so who cares? */
624 greg 2.1 rmx_free(tmp); /* clean up */
625 greg 2.2 rmx_reset(res->rmp);
626     rmx_reset(&res->imx);
627     for (i = 0; i < nmats; i++) {
628     rmx_reset(mop[i].rmp);
629     rmx_reset(&mop[i].imx);
630     if (mop[i].inspec[0] == '!')
631     pclose(mop[i].infp);
632     else if (mop[i].inspec != stdin_name)
633     fclose(mop[i].infp);
634     mop[i].infp = NULL;
635 greg 2.1 }
636 greg 2.2 #endif
637     return(fflush(fout) != EOF);
638 greg 2.1 memerror:
639     fputs("Out of buffer space in combine_input()\n", stderr);
640     return(0);
641     multerror:
642     fputs("Unexpected matrix multiply error in combine_input()\n", stderr);
643     return(0);
644     }
645    
646     static int
647     get_factors(double da[], int n, char *av[])
648     {
649     int ac;
650    
651     for (ac = 0; ac < n && isflt(av[ac]); ac++)
652     da[ac] = atof(av[ac]);
653     return(ac);
654     }
655    
656     static void
657     resize_inparr(int n2alloc)
658     {
659     int i;
660    
661     for (i = nmats; i > n2alloc; i--) {
662     rmx_reset(&mop[i].imx);
663     if (mop[i].rmp != &mop[i].imx)
664     rmx_free(mop[i].rmp);
665     }
666     mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
667     if (mop == NULL) {
668     fputs("Out of memory in resize_inparr()\n", stderr);
669     exit(1);
670     }
671     if (n2alloc > nmats)
672     memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
673     nall = n2alloc;
674     }
675    
676     /* Load one or more matrices and operate on them, sending results to stdout */
677     int
678     main(int argc, char *argv[])
679     {
680    
681     int outfmt = DTfromHeader;
682     const char *defCsym = NULL;
683     int echoheader = 1;
684     int stdin_used = 0;
685     const char *mcat_spec = NULL;
686     int n2comp = 0;
687     uby8 comp_ndx[128];
688     int i;
689     /* get starting input array */
690     mop = (ROPMAT *)calloc(nall=2, sizeof(ROPMAT));
691     /* get options and arguments */
692     for (i = 1; i < argc; i++)
693     if (argv[i][0] != '-' || !argv[i][1]) {
694     if (argv[i][0] == '-') {
695     if (stdin_used++) goto stdin_error;
696     mop[nmats].inspec = stdin_name;
697     } else
698     mop[nmats].inspec = argv[i];
699     if (!mop[nmats].preop.csym)
700     mop[nmats].preop.csym = defCsym;
701     if (++nmats >= nall)
702     resize_inparr(nmats + (nmats>>2) + 2);
703     } else {
704     int n = argc-1 - i;
705     switch (argv[i][1]) { /* get option */
706     case 'w':
707     nowarn = !nowarn;
708     break;
709     case 'h':
710     echoheader = !echoheader;
711     break;
712     case 'e':
713     if (!n) goto userr;
714     comp_ndx[n2comp++] = i++;
715     break;
716     case 'f':
717     switch (argv[i][2]) {
718     case '\0':
719     if (!n) goto userr;
720     comp_ndx[n2comp++] = i++;
721     break;
722     case 'd':
723     outfmt = DTdouble;
724     break;
725     case 'f':
726     outfmt = DTfloat;
727     break;
728     case 'a':
729     outfmt = DTascii;
730     break;
731     case 'c':
732     outfmt = DTrgbe;
733     break;
734     default:
735     goto userr;
736     }
737     break;
738     case 's':
739     if (n > MAXCOMP) n = MAXCOMP;
740     i += mop[nmats].preop.nsf =
741     get_factors(mop[nmats].preop.sca,
742     n, argv+i+1);
743     if (mop[nmats].preop.nsf <= 0) {
744     fprintf(stderr, "%s: -s missing arguments\n",
745     argv[0]);
746     goto userr;
747     }
748     break;
749     case 'C':
750     if (!n || isflt(argv[i+1]))
751     goto userr;
752     defCsym = mop[nmats].preop.csym = argv[++i];
753     mop[nmats].preop.clen = 0;
754     mcat_last = 0;
755     break;
756     case 'c':
757     if (n && !isflt(argv[i+1])) {
758     mop[nmats].preop.csym = argv[++i];
759     mop[nmats].preop.clen = 0;
760     break;
761     }
762     if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
763     i += mop[nmats].preop.clen =
764     get_factors(mop[nmats].preop.cmat,
765     n, argv+i+1);
766     if (mop[nmats].preop.clen <= 0) {
767     fprintf(stderr, "%s: -c missing arguments\n",
768     argv[0]);
769     goto userr;
770     }
771     mop[nmats].preop.csym = NULL;
772     mcat_last = 0;
773     break;
774     case 'm':
775     if (!n) goto userr;
776     if (argv[++i][0] == '-' && !argv[i][1]) {
777     if (stdin_used++) goto stdin_error;
778     mcat_spec = stdin_name;
779     } else
780     mcat_spec = argv[i];
781     mcat_last = 1;
782     break;
783     default:
784     fprintf(stderr, "%s: unknown option '%s'\n",
785     argv[0], argv[i]);
786     goto userr;
787     }
788     }
789     if (!nmats) {
790     fprintf(stderr, "%s: need at least one input matrix\n", argv[0]);
791     goto userr;
792     }
793     resize_inparr(nmats+1); /* extra matrix at end for result */
794     mop[nmats].inspec = "trailing_ops";
795     /* load final concatenation matrix */
796     if (mcat_spec && !(mcat = rmx_load(mcat_spec, RMPnone))) {
797 greg 2.3 fprintf(stderr, "%s: error loading concatenation matrix: %s\n",
798 greg 2.1 argv[0], mcat_spec);
799     return(1);
800     }
801     /* get/check inputs, set constants */
802     if (!initialize(&mop[nmats].imx))
803     return(1);
804    
805     for (i = 0; i < n2comp; i++) /* user .cal files and expressions */
806     if (argv[comp_ndx[i]][1] == 'f') {
807     char *fpath = getpath(argv[comp_ndx[i]+1],
808     getrlibpath(), 0);
809     if (fpath == NULL) {
810     fprintf(stderr, "%s: cannot find file '%s'\n",
811     argv[0], argv[comp_ndx[i]+1]);
812     return(1);
813     }
814     fcompile(fpath);
815     } else /* (argv[comp_ndx[i]][1] == 'e') */
816     scompile(argv[comp_ndx[i]+1], NULL, 0);
817    
818     /* get trailing color transform */
819     if (!get_component_xfm(&mop[nmats]))
820     return(1);
821     /* adjust output dimensions and #components */
822     if (mcat) {
823     if (mop[nmats].imx.ncols != mcat->nrows) {
824     fprintf(stderr,
825     "%s: number of input columns does not match number of rows in '%s'\n",
826     argv[0], mcat_spec);
827     return(1);
828     }
829     if (mcat->ncomp != (mcat_last ? mop[nmats].rmp->ncomp : mop[nmats].imx.ncomp)) {
830     fprintf(stderr,
831     "%s: number of components does not match those in '%s'\n",
832     argv[0], mcat_spec);
833     return(1);
834     }
835 greg 2.3 if (!split_input(&mop[nmats]))
836     return(1);
837 greg 2.1 mop[nmats].rmp->ncols = mcat->ncols;
838     }
839 greg 2.2 newheader("RADIANCE", stdout); /* write output header */
840 greg 2.1 if (echoheader)
841     output_headinfo(stdout);
842     printargs(argc, argv, stdout);
843     fputnow(stdout);
844 greg 2.2 mop[nmats].rmp->dtype = rmx_write_header(mop[nmats].rmp, outfmt, stdout);
845     if (!mop[nmats].rmp->dtype) {
846     fprintf(stderr, "%s: unsupported output format\n", argv[0]);
847 greg 2.1 return(1);
848     }
849     /* process & write rows */
850     return(combine_input(&mop[nmats], stdout) ? 0 : 1);
851     stdin_error:
852     fprintf(stderr, "%s: %s used for more than one input\n",
853     argv[0], stdin_name);
854     return(1);
855     userr:
856     fprintf(stderr,
857     "Usage: %s [-h][-f{adfc}][-e expr][-f file][-s sf .. | -c ce ..] m1 .. -m mcat > mres\n",
858     argv[0]);
859     return(1);
860     }