ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.8
Committed: Thu May 16 18:59:19 2024 UTC (23 hours, 45 minutes ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.7: +1 -2 lines
Log Message:
fix: Made use of resolu.h more consistent and reliable

File Contents

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