ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.7
Committed: Fri Feb 23 03:47:57 2024 UTC (2 months, 3 weeks ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.6: +2 -1 lines
Log Message:
perf: Added array index optimization to calcomp routines

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.7 static const char RCSid[] = "$Id: rcomb.c,v 2.6 2024/02/09 19:58:56 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 greg 2.2 int in_nrows; /* number of input rows (or 0) */
49     #define in_ncols (mop[0].rmp->ncols) /* number of input columns */
50 greg 2.1 #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     if (refm.ncomp == 3) {
93     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 greg 2.3 double cf = 1;
190 greg 2.1 int i, j;
191     /* check suffix => reference file */
192     if (strchr(rop->preop.csym, '.') > rop->preop.csym)
193     return(checkreffile(rop));
194    
195     if (nc < 3) {
196     fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
197     rop->inspec, rop->preop.csym);
198     return(0);
199     }
200     rop->preop.clen = strlen(rop->preop.csym) * nc;
201     if (rop->preop.clen > MAXCOMP*MAXCOMP) {
202     fprintf(stderr, "%s: -c '%s' results in too many components\n",
203     rop->inspec, rop->preop.csym);
204     return(0);
205     }
206     for (j = 0; rop->preop.csym[j]; j++) {
207     int comp = 0;
208     switch (rop->preop.csym[j]) {
209     case 'B':
210 greg 2.3 case 'b':
211 greg 2.1 ++comp;
212     /* fall through */
213     case 'G':
214 greg 2.3 case 'g':
215 greg 2.1 ++comp;
216     /* fall through */
217     case 'R':
218 greg 2.3 case 'r':
219     if (rop->preop.csym[j] <= 'Z')
220     cf = 1./WHTEFFICACY;
221 greg 2.1 if (dt == DTxyze) {
222     for (i = 3; i--; )
223 greg 2.3 rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
224 greg 2.1 } else if (nc == 3)
225     rop->preop.cmat[j*nc+comp] = 1.;
226     else
227     rgbrow(rop, j, comp);
228     break;
229     case 'Z':
230 greg 2.3 case 'z':
231 greg 2.1 ++comp;
232     /* fall through */
233     case 'Y':
234 greg 2.3 case 'y':
235 greg 2.1 ++comp;
236     /* fall through */
237     case 'X':
238 greg 2.3 case 'x':
239     if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
240     cf = WHTEFFICACY;
241 greg 2.1 if (dt == DTxyze) {
242     rop->preop.cmat[j*nc+comp] = 1.;
243     } else if (nc == 3) {
244     for (i = 3; i--; )
245     rop->preop.cmat[j*nc+i] =
246     rgb2xyzmat[comp][i];
247     } else if (comp == CIEY)
248     sensrow(rop, j, scolor2photopic);
249     else
250     xyzrow(rop, j, comp);
251    
252 greg 2.3 for (i = nc*(cf != 1); i--; )
253     rop->preop.cmat[j*nc+i] *= cf;
254 greg 2.1 break;
255     case 'S': /* scotopic (il)luminance */
256 greg 2.3 cf = WHTSCOTOPIC;
257     /* fall through */
258     case 's':
259 greg 2.1 sensrow(rop, j, scolor2scotopic);
260 greg 2.3 for (i = nc*(cf != 1); i--; )
261     rop->preop.cmat[j*nc+i] *= cf;
262 greg 2.1 break;
263     case 'M': /* melanopic (il)luminance */
264 greg 2.3 cf = WHTMELANOPIC;
265     /* fall through */
266     case 'm':
267 greg 2.1 sensrow(rop, j, scolor2melanopic);
268 greg 2.3 for (i = nc*(cf != 1); i--; )
269     rop->preop.cmat[j*nc+i] *= cf;
270 greg 2.1 break;
271     case 'A': /* average component */
272 greg 2.3 case 'a':
273 greg 2.1 for (i = nc; i--; )
274     rop->preop.cmat[j*nc+i] = 1./(double)nc;
275     break;
276     default:
277     fprintf(stderr, "%s: -c '%c' unsupported\n",
278     rop->inspec, rop->preop.csym[j]);
279     return(0);
280     }
281     }
282     if (!split_input(rop)) /* get our own struct */
283     return(0);
284     memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
285     rop->rmp->ncomp = rop->preop.clen / nc;
286     /* decide on output type */
287 greg 2.4 if (!strcasecmp(rop->preop.csym, "XYZ")) {
288 greg 2.1 if (dt <= DTspec)
289     rop->rmp->dtype = DTxyze;
290 greg 2.4 } else if (!strcasecmp(rop->preop.csym, "RGB")) {
291 greg 2.1 if (dt <= DTspec)
292     rop->rmp->dtype = DTrgbe;
293     } else if (rop->rmp->dtype == DTspec)
294     rop->rmp->dtype = DTfloat;
295     return(1);
296     }
297    
298     static int
299     get_component_xfm(ROPMAT *rop)
300     {
301     int i, j;
302    
303     if (rop->rmp != &rop->imx) { /* reset destination matrix */
304     rmx_free(rop->rmp);
305     rop->rmp = &rop->imx;
306     }
307     if (rop->preop.csym && /* symbolic transform? */
308     !checksymbolic(rop))
309     return(0);
310     /* undo exposure? */
311     if (fabs(1. - bright(rop->rmp->cexp)) > .025) {
312     if (rop->rmp->ncomp == 1)
313     rop->rmp->cexp[RED] = rop->rmp->cexp[GRN] =
314     rop->rmp->cexp[BLU] = bright(rop->rmp->cexp);
315     if (rop->preop.nsf <= 0) {
316     rop->preop.nsf = i = rop->rmp->ncomp;
317     while (i--)
318     rop->preop.sca[i] = 1.;
319     }
320     if (rop->preop.nsf == 1) {
321     if (rop->rmp->ncomp == 3) {
322     rop->preop.sca[2] = rop->preop.sca[1] =
323     rop->preop.sca[0];
324     rop->preop.nsf = 3;
325     } else
326     rop->preop.sca[0] /= bright(rop->rmp->cexp);
327     }
328     if (rop->preop.nsf == 3) {
329     opcolor(rop->preop.sca, /=, rop->rmp->cexp);
330     } else if (rop->preop.nsf > 3) { /* punt */
331     double mult = 1./bright(rop->rmp->cexp);
332     for (i = rop->preop.nsf; i--; )
333     rop->preop.sca[i] *= mult;
334     }
335     setcolor(rop->rmp->cexp, 1., 1., 1.);
336     }
337     if (rop->preop.clen > 0) { /* use component transform? */
338     if (rop->preop.clen % rop->imx.ncomp) {
339     fprintf(stderr, "%s: -c must have N x %d coefficients\n",
340     rop->inspec, rop->imx.ncomp);
341     return(0);
342     }
343     if (rop->preop.nsf > 0) { /* scale transform, instead */
344     if (rop->preop.nsf == 1) {
345     for (i = rop->preop.clen; i--; )
346     rop->preop.cmat[i] *= rop->preop.sca[0];
347     } else if (rop->preop.nsf*rop->imx.ncomp != rop->preop.clen) {
348     fprintf(stderr, "%s: -s must have one or %d factors\n",
349     rop->inspec,
350     rop->preop.clen/rop->imx.ncomp);
351     return(0);
352     } else {
353     for (i = rop->preop.nsf; i--; )
354     for (j = rop->imx.ncomp; j--; )
355     rop->preop.cmat[i*rop->imx.ncomp+j]
356     *= rop->preop.sca[i];
357     }
358     }
359     rop->preop.nsf = 0; /* now folded in */
360     if (!split_input(rop)) /* get our own struct */
361     return(0);
362     rop->rmp->ncomp = rop->preop.clen / rop->imx.ncomp;
363     if ((rop->rmp->ncomp > 3) & (rop->rmp->dtype <= DTspec)) {
364     rop->rmp->dtype = DTfloat; /* probably not actual spectrum */
365     memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
366     }
367     } else if (rop->preop.nsf > 0) { /* else use scalar(s)? */
368     if (rop->preop.nsf == 1) {
369     for (i = rop->rmp->ncomp; --i; )
370     rop->preop.sca[i] = rop->preop.sca[0];
371     rop->preop.nsf = rop->rmp->ncomp;
372     } else if (rop->preop.nsf != rop->rmp->ncomp) {
373     fprintf(stderr, "%s: -s must have one or %d factors\n",
374     rop->inspec, rop->rmp->ncomp);
375     return(0);
376     }
377     }
378     return(1);
379     }
380    
381     static int
382     apply_op(RMATRIX *dst, const RMATRIX *src, const RUNARYOP *ro)
383     {
384     if (ro->clen > 0) {
385     RMATRIX *res = rmx_transform(src, dst->ncomp, ro->cmat);
386     if (!res) {
387     fputs("Error in call to rmx_transform()\n", stderr);
388     return(0);
389     }
390     if (!rmx_transfer_data(dst, res, 0))
391     return(0);
392     rmx_free(res);
393     } else if (dst != src)
394     memcpy(dst->mtx, src->mtx,
395     sizeof(double)*dst->ncomp*dst->ncols*dst->nrows);
396     if (ro->nsf == dst->ncomp)
397     rmx_scale(dst, ro->sca);
398     return(1);
399     }
400    
401     static int
402     open_input(ROPMAT *rop)
403     {
404     int outtype;
405    
406     if (!rop || !rop->inspec || !rop->inspec[0])
407     return(0);
408     if (rop->inspec == stdin_name)
409     rop->infp = stdin;
410     else if (rop->inspec[0] == '!')
411     rop->infp = popen(rop->inspec+1, "r");
412     else
413     rop->infp = fopen(rop->inspec, "rb");
414    
415     if (!rmx_load_header(&rop->imx, rop->infp)) {
416     fprintf(stderr, "Bad header from: %s\n", rop->inspec);
417     return(0);
418     }
419     return(get_component_xfm(rop));
420     }
421    
422     /* Return nominal wavelength associated with input component (return nm) */
423     static double
424     l_wavelength(char *nam)
425     {
426     double comp = argument(1);
427    
428     if ((comp < -.5) | (comp >= in_ncomp+.5)) {
429     errno = EDOM;
430     return(.0);
431     }
432     if (comp < .5) /* asking for #components? */
433     return(in_ncomp);
434    
435     if (in_ncomp == 3) { /* special case for RGB */
436     const int w0 = (int)(comp - .5);
437     return(mop[0].rmp->wlpart[w0] +
438     (comp-.5)*(mop[0].rmp->wlpart[w0+1] -
439     mop[0].rmp->wlpart[w0]));
440     }
441     return(mop[0].rmp->wlpart[0] + /* general case, even div. */
442     (comp-.5)/(double)in_ncomp *
443     (mop[0].rmp->wlpart[3] - mop[0].rmp->wlpart[0]));
444     }
445    
446     /* Return ith input with optional channel selector */
447     static double
448     l_chanin(char *nam)
449     {
450     double inp = argument(1);
451     int mi, chan;
452    
453     if ((mi = (int)(inp-.5)) < 0 || mi >= nmats) {
454     errno = EDOM;
455     return(.0);
456     }
457     if (inp < .5) /* asking for #inputs? */
458     return(nmats);
459    
460     if (nargum() >= 2) {
461     double cval = argument(2);
462     if (cval < .5 || (chan = (int)(cval-.5)) >= in_ncomp) {
463     errno = EDOM;
464     return(.0);
465     }
466     } else
467     chan = cur_chan;
468    
469     return(mop[mi].rmp->mtx[cur_col*in_ncomp + chan]);
470     }
471    
472     static int
473     initialize(RMATRIX *imp)
474     {
475     int i;
476     /* XXX struct is zeroed coming in */
477     setcolor(imp->cexp, 1.f, 1.f, 1.f);
478     for (i = 0; i < nmats; i++) { /* open each input */
479     int restype;
480     if (!open_input(&mop[i]))
481     return(0);
482     restype = mop[i].rmp->dtype;
483     if (!imp->dtype || (restype = rmx_newtype(restype, imp->dtype)) > 0)
484     imp->dtype = restype;
485     else
486     fprintf(stderr, "%s: warning - data type mismatch\n",
487     mop[i].inspec);
488     if (!i) {
489     imp->ncols = mop[0].rmp->ncols;
490     imp->ncomp = mop[0].rmp->ncomp;
491     memcpy(imp->wlpart, mop[0].rmp->wlpart, sizeof(imp->wlpart));
492 greg 2.2 } else if ((mop[i].rmp->ncols != imp->ncols) |
493     (mop[i].rmp->ncomp != imp->ncomp) |
494     ((in_nrows > 0) & (mop[i].rmp->nrows > 0) &
495     (mop[i].rmp->nrows != in_nrows))) {
496 greg 2.1 fprintf(stderr, "%s: mismatch in size or #components\n",
497     mop[i].inspec);
498     return(0);
499     } /* XXX should check wlpart? */
500 greg 2.2 if (in_nrows <= 0)
501     in_nrows = imp->nrows = mop[i].rmp->nrows;
502 greg 2.1 } /* set up .cal environment */
503     esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
504     esupport &= ~(E_OUTCHAN|E_INCHAN);
505     varset("PI", ':', M_PI);
506     varset("nfiles", ':', nmats);
507     varset("nrows", ':', in_nrows);
508     varset("ncols", ':', in_ncols);
509     varset("ncomp", ':', in_ncomp);
510     varset("R", ':', 1.);
511     varset("G", ':', 2.);
512     varset("B", ':', 3.);
513     funset("wl", 1, ':', l_wavelength);
514     funset("ci", 1, '=', l_chanin);
515     scompile("ri(i)=ci(i,R);gi(i)=ci(i,G);bi(i)=ci(i,B)", NULL, 0);
516     return(1);
517     }
518    
519     static void
520     output_headinfo(FILE *fp)
521     {
522     int i;
523    
524     for (i = 0; i < nmats; i++) {
525     const char *cp = mop[i].imx.info;
526     fputs(mop[i].inspec, fp);
527     fputs(":\n", fp);
528     if (!cp) continue;
529     while (*cp) {
530     if (*cp == '\n') {
531     cp++; /* avoid inadvertant terminus */
532     continue;
533     }
534     fputc('\t', fp); /* indent this input's info */
535     do
536     putc(*cp, fp);
537     while (*cp++ != '\n');
538     }
539     }
540     }
541    
542     static int
543     combine_input(ROPMAT *res, FILE *fout)
544     {
545     int set_r, set_c;
546     RMATRIX *tmp = NULL;
547     int co_set;
548     int i;
549     /* allocate input row buffers */
550     for (i = 0; i < nmats; i++) {
551     mop[i].imx.nrows = 1; /* we'll be doing a row at a time */
552     if (!rmx_prepare(&mop[i].imx))
553     goto memerror;
554     if (mop[i].rmp != &mop[i].imx) {
555     mop[i].rmp->nrows = 1;
556     if (!rmx_prepare(mop[i].rmp))
557     goto memerror;
558     }
559     }
560     /* prep output row buffers */
561     if (mcat || res->preop.clen > 0) {
562     if (!split_input(res)) /* need separate buffer */
563     return(0);
564     if (res->preop.clen > 0)
565     res->rmp->ncomp = res->preop.clen / res->imx.ncomp;
566     res->rmp->nrows = 1;
567     if (!mcat | !mcat_last && !rmx_prepare(res->rmp))
568     goto memerror;
569     }
570     if (mcat && mcat_last &&
571     !(tmp = rmx_alloc(1, res->imx.ncols, res->rmp->ncomp)))
572     goto memerror;
573     res->imx.nrows = 1;
574     if (!rmx_prepare(&res->imx))
575     goto memerror;
576     /* figure out what the user set */
577     co_set = fundefined("co");
578     if (!co_set)
579     co_set = -vardefined("co");
580     if (!co_set & (in_ncomp == 3) && vardefined("ro") &&
581     vardefined("go") && vardefined("bo")) {
582     scompile("co(p)=select(p,ro,go,bo)", NULL, 0);
583     co_set = 1;
584     }
585     if (co_set) { /* set if user wants, didn't set */
586     set_r = varlookup("r") != NULL && !vardefined("r");
587     set_c = varlookup("c") != NULL && !vardefined("c");
588     } else /* save a little time */
589     set_r = set_c = 0;
590     /* read/process row-by-row */
591 greg 2.2 for (cur_row = 0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row++) {
592 greg 2.1 RMATRIX *mres = NULL;
593     for (i = 0; i < nmats; i++) {
594     if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
595 greg 2.5 if (cur_row > in_nrows) /* unknown #input rows? */
596 greg 2.2 goto loop_exit;
597 greg 2.1 fprintf(stderr, "%s: read error at row %d\n",
598     mop[i].inspec, cur_row);
599     return(0);
600     }
601     if (!apply_op(mop[i].rmp, &mop[i].imx, &mop[i].preop))
602     return(0);
603     }
604     if (set_r) varset("r", '=', cur_row);
605     for (cur_col = 0; cur_col < in_ncols; cur_col++) {
606     if (set_c) varset("c", '=', cur_col);
607     for (cur_chan = 0; cur_chan < in_ncomp; cur_chan++) {
608     const int ndx = cur_col*in_ncomp + cur_chan;
609     eclock++;
610     if (!co_set) { /* just summing elements? */
611     res->imx.mtx[ndx] = 0;
612     for (i = nmats; i--; )
613     res->imx.mtx[ndx] += mop[i].rmp->mtx[ndx];
614     } else if (co_set > 0) {
615     double dchan = cur_chan+1;
616     res->imx.mtx[ndx] = funvalue("co", 1, &dchan);
617     } else
618     res->imx.mtx[ndx] = varvalue("co");
619     }
620     } /* final conversions */
621     if (!mcat) {
622     if (!apply_op(res->rmp, &res->imx, &res->preop))
623     return(0);
624     } else if (mcat_last) {
625     if (!apply_op(tmp, &res->imx, &res->preop))
626     return(0);
627     mres = rmx_multiply(tmp, mcat);
628     if (!mres)
629     goto multerror;
630     if (!rmx_transfer_data(res->rmp, mres, 0))
631     return(0);
632     } else /* mcat && !mcat_last */ {
633     mres = rmx_multiply(&res->imx, mcat);
634     if (!mres)
635     goto multerror;
636     if (!apply_op(res->rmp, mres, &res->preop))
637     return(0);
638     }
639     rmx_free(mres); mres = NULL;
640     if (!rmx_write_data(res->rmp->mtx, res->rmp->ncomp,
641     res->rmp->ncols, res->rmp->dtype, fout))
642     return(0);
643     }
644 greg 2.2 loop_exit:
645 greg 2.1 #if 0 /* we're about to exit, so who cares? */
646     rmx_free(tmp); /* clean up */
647     rmx_reset(res->rmp);
648     rmx_reset(&res->imx);
649     for (i = 0; i < nmats; i++) {
650     rmx_reset(mop[i].rmp);
651     rmx_reset(&mop[i].imx);
652     if (mop[i].inspec[0] == '!')
653     pclose(mop[i].infp);
654     else if (mop[i].inspec != stdin_name)
655     fclose(mop[i].infp);
656     mop[i].infp = NULL;
657     }
658     #endif
659     return(fflush(fout) != EOF);
660     memerror:
661     fputs("Out of buffer space in combine_input()\n", stderr);
662     return(0);
663     multerror:
664     fputs("Unexpected matrix multiply error in combine_input()\n", stderr);
665     return(0);
666     }
667    
668     static int
669     get_factors(double da[], int n, char *av[])
670     {
671     int ac;
672    
673     for (ac = 0; ac < n && isflt(av[ac]); ac++)
674     da[ac] = atof(av[ac]);
675     return(ac);
676     }
677    
678     static void
679     resize_inparr(int n2alloc)
680     {
681     int i;
682    
683 greg 2.6 if (n2alloc == nall)
684     return;
685     for (i = nall; i > n2alloc; i--) {
686 greg 2.1 rmx_reset(&mop[i].imx);
687     if (mop[i].rmp != &mop[i].imx)
688     rmx_free(mop[i].rmp);
689     }
690     mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
691     if (mop == NULL) {
692     fputs("Out of memory in resize_inparr()\n", stderr);
693     exit(1);
694     }
695 greg 2.6 if (n2alloc > nall)
696     memset(mop+nall, 0, (n2alloc-nall)*sizeof(ROPMAT));
697 greg 2.1 nall = n2alloc;
698     }
699    
700     /* Load one or more matrices and operate on them, sending results to stdout */
701     int
702     main(int argc, char *argv[])
703     {
704    
705     int outfmt = DTfromHeader;
706     const char *defCsym = NULL;
707     int echoheader = 1;
708     int stdin_used = 0;
709     const char *mcat_spec = NULL;
710     int n2comp = 0;
711     uby8 comp_ndx[128];
712     int i;
713     /* get starting input array */
714     mop = (ROPMAT *)calloc(nall=2, sizeof(ROPMAT));
715     /* get options and arguments */
716     for (i = 1; i < argc; i++)
717     if (argv[i][0] != '-' || !argv[i][1]) {
718     if (argv[i][0] == '-') {
719     if (stdin_used++) goto stdin_error;
720     mop[nmats].inspec = stdin_name;
721     } else
722     mop[nmats].inspec = argv[i];
723     if (!mop[nmats].preop.csym)
724     mop[nmats].preop.csym = defCsym;
725     if (++nmats >= nall)
726     resize_inparr(nmats + (nmats>>2) + 2);
727     } else {
728     int n = argc-1 - i;
729     switch (argv[i][1]) { /* get option */
730     case 'w':
731     nowarn = !nowarn;
732     break;
733     case 'h':
734     echoheader = !echoheader;
735     break;
736     case 'e':
737     if (!n) goto userr;
738     comp_ndx[n2comp++] = i++;
739     break;
740     case 'f':
741     switch (argv[i][2]) {
742     case '\0':
743     if (!n) goto userr;
744     comp_ndx[n2comp++] = i++;
745     break;
746     case 'd':
747     outfmt = DTdouble;
748     break;
749     case 'f':
750     outfmt = DTfloat;
751     break;
752     case 'a':
753     outfmt = DTascii;
754     break;
755     case 'c':
756     outfmt = DTrgbe;
757     break;
758     default:
759     goto userr;
760     }
761     break;
762     case 's':
763     if (n > MAXCOMP) n = MAXCOMP;
764     i += mop[nmats].preop.nsf =
765     get_factors(mop[nmats].preop.sca,
766     n, argv+i+1);
767     if (mop[nmats].preop.nsf <= 0) {
768     fprintf(stderr, "%s: -s missing arguments\n",
769     argv[0]);
770     goto userr;
771     }
772     break;
773     case 'C':
774     if (!n || isflt(argv[i+1]))
775     goto userr;
776     defCsym = mop[nmats].preop.csym = argv[++i];
777     mop[nmats].preop.clen = 0;
778     mcat_last = 0;
779     break;
780     case 'c':
781     if (n && !isflt(argv[i+1])) {
782     mop[nmats].preop.csym = argv[++i];
783     mop[nmats].preop.clen = 0;
784     break;
785     }
786     if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
787     i += mop[nmats].preop.clen =
788     get_factors(mop[nmats].preop.cmat,
789     n, argv+i+1);
790     if (mop[nmats].preop.clen <= 0) {
791     fprintf(stderr, "%s: -c missing arguments\n",
792     argv[0]);
793     goto userr;
794     }
795     mop[nmats].preop.csym = NULL;
796     mcat_last = 0;
797     break;
798     case 'm':
799     if (!n) goto userr;
800     if (argv[++i][0] == '-' && !argv[i][1]) {
801     if (stdin_used++) goto stdin_error;
802     mcat_spec = stdin_name;
803     } else
804     mcat_spec = argv[i];
805     mcat_last = 1;
806     break;
807     default:
808     fprintf(stderr, "%s: unknown option '%s'\n",
809     argv[0], argv[i]);
810     goto userr;
811     }
812     }
813     if (!nmats) {
814     fprintf(stderr, "%s: need at least one input matrix\n", argv[0]);
815     goto userr;
816     }
817     resize_inparr(nmats+1); /* extra matrix at end for result */
818     mop[nmats].inspec = "trailing_ops";
819     /* load final concatenation matrix */
820     if (mcat_spec && !(mcat = rmx_load(mcat_spec, RMPnone))) {
821     fprintf(stderr, "%s: error loading concatenation matrix: %s\n",
822     argv[0], mcat_spec);
823     return(1);
824     }
825     /* get/check inputs, set constants */
826     if (!initialize(&mop[nmats].imx))
827     return(1);
828    
829     for (i = 0; i < n2comp; i++) /* user .cal files and expressions */
830     if (argv[comp_ndx[i]][1] == 'f') {
831     char *fpath = getpath(argv[comp_ndx[i]+1],
832     getrlibpath(), 0);
833     if (fpath == NULL) {
834     fprintf(stderr, "%s: cannot find file '%s'\n",
835     argv[0], argv[comp_ndx[i]+1]);
836     return(1);
837     }
838     fcompile(fpath);
839     } else /* (argv[comp_ndx[i]][1] == 'e') */
840     scompile(argv[comp_ndx[i]+1], NULL, 0);
841    
842     /* get trailing color transform */
843     if (!get_component_xfm(&mop[nmats]))
844     return(1);
845     /* adjust output dimensions and #components */
846     if (mcat) {
847     if (mop[nmats].imx.ncols != mcat->nrows) {
848     fprintf(stderr,
849     "%s: number of input columns does not match number of rows in '%s'\n",
850     argv[0], mcat_spec);
851     return(1);
852     }
853     if (mcat->ncomp != (mcat_last ? mop[nmats].rmp->ncomp : mop[nmats].imx.ncomp)) {
854     fprintf(stderr,
855     "%s: number of components does not match those in '%s'\n",
856     argv[0], mcat_spec);
857     return(1);
858     }
859     if (!split_input(&mop[nmats]))
860     return(1);
861     mop[nmats].rmp->ncols = mcat->ncols;
862     }
863     newheader("RADIANCE", stdout); /* write output header */
864     if (echoheader)
865     output_headinfo(stdout);
866     printargs(argc, argv, stdout);
867     fputnow(stdout);
868     mop[nmats].rmp->dtype = rmx_write_header(mop[nmats].rmp, outfmt, stdout);
869     if (!mop[nmats].rmp->dtype) {
870     fprintf(stderr, "%s: unsupported output format\n", argv[0]);
871     return(1);
872     }
873 greg 2.7 doptimize(1); /* optimize definitions */
874 greg 2.1 /* process & write rows */
875     return(combine_input(&mop[nmats], stdout) ? 0 : 1);
876     stdin_error:
877     fprintf(stderr, "%s: %s used for more than one input\n",
878     argv[0], stdin_name);
879     return(1);
880     userr:
881     fprintf(stderr,
882     "Usage: %s [-h][-f{adfc}][-e expr][-f file][-s sf .. | -c ce ..] m1 .. -m mcat > mres\n",
883     argv[0]);
884     return(1);
885     }