ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.15
Committed: Thu May 23 01:28:07 2024 UTC (3 weeks, 3 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.14: +2 -1 lines
Log Message:
fix(rcomb): Fixed confusion over #components with symbolic -c in final xform

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.15 static const char RCSid[] = "$Id: rcomb.c,v 2.14 2024/05/22 15:38:04 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * General component matrix combiner, operating on a row at a time.
6     */
7    
8 greg 2.12 #include <signal.h>
9 greg 2.1 #include <math.h>
10     #include "platform.h"
11 greg 2.10 #include "rtprocess.h"
12 greg 2.1 #include "rtio.h"
13     #include "rmatrix.h"
14     #include "calcomp.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 greg 2.10 SUBPROC *cproc = NULL; /* child process array */
58     int nchildren = 0; /* # of child processes */
59     int inchild = -1; /* our child ID (-1: parent) */
60 greg 2.13 int pgid = -1; /* process group ID */
61 greg 2.12 int nr_out = 0; /* # of rows output by kids */
62 greg 2.10
63 greg 2.1 static int checksymbolic(ROPMAT *rop);
64    
65 greg 2.12 static void
66     on_sigio(int dummy)
67     {
68     nr_out++; /* happens when child outputs row */
69     }
70    
71 greg 2.1 static int
72     split_input(ROPMAT *rop)
73     {
74     if (rop->rmp == &rop->imx && !(rop->rmp = rmx_copy(&rop->imx))) {
75     fputs("Out of memory in split_input()\n", stderr);
76     return(0);
77     }
78     rmx_reset(rop->rmp);
79     return(1);
80     }
81    
82     /* Check/set transform based on a reference input file */
83     static int
84     checkreffile(ROPMAT *rop)
85     {
86     static const char *curRF = NULL;
87     static RMATRIX refm;
88     const int nc = rop->imx.ncomp;
89     int i;
90    
91     if (!curRF || strcmp(rop->preop.csym, curRF)) {
92     FILE *fp = fopen(rop->preop.csym, "rb");
93     if (!rmx_load_header(&refm, fp)) {
94     fprintf(stderr, "%s: cannot read info header\n",
95     rop->preop.csym);
96     curRF = NULL;
97     if (fp) fclose(fp);
98     return(0);
99     }
100     fclose(fp);
101     curRF = rop->preop.csym;
102     }
103     if (refm.ncomp == 3) {
104     rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
105     return(checksymbolic(rop));
106     }
107     if (refm.ncomp == 2) {
108     fprintf(stderr, "%s: cannot convert to 2 components\n",
109     curRF);
110     return(0);
111     }
112     if (refm.ncomp == 1) {
113     rop->preop.csym = "Y"; /* XXX big assumption */
114     return(checksymbolic(rop));
115     }
116     if (refm.ncomp == nc &&
117     !memcmp(refm.wlpart, rop->imx.wlpart, sizeof(refm.wlpart)))
118     return(1); /* nothing to do */
119    
120     if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
121     fprintf(stderr, "%s: cannot resample from %d to %d components\n",
122     curRF, nc, refm.ncomp);
123     return(0);
124     }
125     if (!split_input(rop)) /* get our own struct */
126     return(0);
127     rop->preop.clen = refm.ncomp * nc; /* compute spec to ref */
128    
129     for (i = 0; i < nc; i++) {
130     SCOLOR scstim, scresp;
131     int j;
132     memset(scstim, 0, sizeof(COLORV)*nc);
133     scstim[i] = 1.f;
134     convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
135     scstim, nc, rop->imx.wlpart[0], rop->imx.wlpart[3]);
136     for (j = refm.ncomp; j-- > 0; )
137     rop->preop.cmat[j*nc + i] = scresp[j];
138     }
139     /* remember new spectral params */
140     memcpy(rop->rmp->wlpart, refm.wlpart, sizeof(rop->rmp->wlpart));
141     rop->rmp->ncomp = refm.ncomp;
142     return(1);
143     }
144    
145     /* Compute conversion row from spectrum to one channel of RGB */
146     static void
147     rgbrow(ROPMAT *rop, int r, int p)
148     {
149     const int nc = rop->imx.ncomp;
150     const float * wlp = rop->imx.wlpart;
151     int i;
152    
153     for (i = nc; i--; ) {
154     int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
155     int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
156     COLOR crgb;
157     spec_rgb(crgb, nmStart, nmEnd);
158     rop->preop.cmat[r*nc+i] = crgb[p];
159     }
160     }
161    
162     /* Compute conversion row from spectrum to one channel of XYZ */
163     static void
164     xyzrow(ROPMAT *rop, int r, int p)
165     {
166     const int nc = rop->imx.ncomp;
167     const float * wlp = rop->imx.wlpart;
168     int i;
169    
170     for (i = nc; i--; ) {
171     int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
172     int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
173     COLOR cxyz;
174     spec_cie(cxyz, nmStart, nmEnd);
175     rop->preop.cmat[r*nc+i] = cxyz[p];
176     }
177     }
178    
179     /* Use the spectral sensitivity function to compute matrix coefficients */
180     static void
181     sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
182     {
183     const int nc = rop->imx.ncomp;
184     int i;
185    
186     for (i = nc; i--; ) {
187     SCOLOR sclr;
188     memset(sclr, 0, sizeof(COLORV)*nc);
189     sclr[i] = 1.f;
190     rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->imx.wlpart);
191     }
192     }
193    
194     /* Check/set symbolic transform */
195     static int
196     checksymbolic(ROPMAT *rop)
197     {
198     const int nc = rop->imx.ncomp;
199     const int dt = rop->imx.dtype;
200 greg 2.3 double cf = 1;
201 greg 2.1 int i, j;
202     /* check suffix => reference file */
203     if (strchr(rop->preop.csym, '.') > rop->preop.csym)
204     return(checkreffile(rop));
205    
206     if (nc < 3) {
207     fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
208     rop->inspec, rop->preop.csym);
209     return(0);
210     }
211     rop->preop.clen = strlen(rop->preop.csym) * nc;
212     if (rop->preop.clen > MAXCOMP*MAXCOMP) {
213     fprintf(stderr, "%s: -c '%s' results in too many components\n",
214     rop->inspec, rop->preop.csym);
215     return(0);
216     }
217     for (j = 0; rop->preop.csym[j]; j++) {
218     int comp = 0;
219     switch (rop->preop.csym[j]) {
220     case 'B':
221 greg 2.3 case 'b':
222 greg 2.1 ++comp;
223     /* fall through */
224     case 'G':
225 greg 2.3 case 'g':
226 greg 2.1 ++comp;
227     /* fall through */
228     case 'R':
229 greg 2.3 case 'r':
230     if (rop->preop.csym[j] <= 'Z')
231     cf = 1./WHTEFFICACY;
232 greg 2.1 if (dt == DTxyze) {
233     for (i = 3; i--; )
234 greg 2.3 rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
235 greg 2.1 } else if (nc == 3)
236     rop->preop.cmat[j*nc+comp] = 1.;
237     else
238     rgbrow(rop, j, comp);
239     break;
240     case 'Z':
241 greg 2.3 case 'z':
242 greg 2.1 ++comp;
243     /* fall through */
244     case 'Y':
245 greg 2.3 case 'y':
246 greg 2.1 ++comp;
247     /* fall through */
248     case 'X':
249 greg 2.3 case 'x':
250     if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
251     cf = WHTEFFICACY;
252 greg 2.1 if (dt == DTxyze) {
253     rop->preop.cmat[j*nc+comp] = 1.;
254     } else if (nc == 3) {
255     for (i = 3; i--; )
256     rop->preop.cmat[j*nc+i] =
257     rgb2xyzmat[comp][i];
258     } else if (comp == CIEY)
259     sensrow(rop, j, scolor2photopic);
260     else
261     xyzrow(rop, j, comp);
262    
263 greg 2.3 for (i = nc*(cf != 1); i--; )
264     rop->preop.cmat[j*nc+i] *= cf;
265 greg 2.1 break;
266     case 'S': /* scotopic (il)luminance */
267 greg 2.3 cf = WHTSCOTOPIC;
268     /* fall through */
269     case 's':
270 greg 2.1 sensrow(rop, j, scolor2scotopic);
271 greg 2.3 for (i = nc*(cf != 1); i--; )
272     rop->preop.cmat[j*nc+i] *= cf;
273 greg 2.1 break;
274     case 'M': /* melanopic (il)luminance */
275 greg 2.3 cf = WHTMELANOPIC;
276     /* fall through */
277     case 'm':
278 greg 2.1 sensrow(rop, j, scolor2melanopic);
279 greg 2.3 for (i = nc*(cf != 1); i--; )
280     rop->preop.cmat[j*nc+i] *= cf;
281 greg 2.1 break;
282     case 'A': /* average component */
283 greg 2.3 case 'a':
284 greg 2.1 for (i = nc; i--; )
285     rop->preop.cmat[j*nc+i] = 1./(double)nc;
286     break;
287     default:
288     fprintf(stderr, "%s: -c '%c' unsupported\n",
289     rop->inspec, rop->preop.csym[j]);
290     return(0);
291     }
292     }
293     if (!split_input(rop)) /* get our own struct */
294     return(0);
295     memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
296     rop->rmp->ncomp = rop->preop.clen / nc;
297     /* decide on output type */
298 greg 2.4 if (!strcasecmp(rop->preop.csym, "XYZ")) {
299 greg 2.1 if (dt <= DTspec)
300     rop->rmp->dtype = DTxyze;
301 greg 2.4 } else if (!strcasecmp(rop->preop.csym, "RGB")) {
302 greg 2.1 if (dt <= DTspec)
303     rop->rmp->dtype = DTrgbe;
304     } else if (rop->rmp->dtype == DTspec)
305     rop->rmp->dtype = DTfloat;
306     return(1);
307     }
308    
309     static int
310     get_component_xfm(ROPMAT *rop)
311     {
312     int i, j;
313    
314     if (rop->rmp != &rop->imx) { /* reset destination matrix */
315     rmx_free(rop->rmp);
316     rop->rmp = &rop->imx;
317     }
318     if (rop->preop.csym && /* symbolic transform? */
319     !checksymbolic(rop))
320     return(0);
321     /* undo exposure? */
322     if (fabs(1. - bright(rop->rmp->cexp)) > .025) {
323     if (rop->rmp->ncomp == 1)
324     rop->rmp->cexp[RED] = rop->rmp->cexp[GRN] =
325     rop->rmp->cexp[BLU] = bright(rop->rmp->cexp);
326     if (rop->preop.nsf <= 0) {
327     rop->preop.nsf = i = rop->rmp->ncomp;
328     while (i--)
329     rop->preop.sca[i] = 1.;
330     }
331     if (rop->preop.nsf == 1) {
332     if (rop->rmp->ncomp == 3) {
333     rop->preop.sca[2] = rop->preop.sca[1] =
334     rop->preop.sca[0];
335     rop->preop.nsf = 3;
336     } else
337     rop->preop.sca[0] /= bright(rop->rmp->cexp);
338     }
339     if (rop->preop.nsf == 3) {
340     opcolor(rop->preop.sca, /=, rop->rmp->cexp);
341     } else if (rop->preop.nsf > 3) { /* punt */
342     double mult = 1./bright(rop->rmp->cexp);
343     for (i = rop->preop.nsf; i--; )
344     rop->preop.sca[i] *= mult;
345     }
346     setcolor(rop->rmp->cexp, 1., 1., 1.);
347     }
348     if (rop->preop.clen > 0) { /* use component transform? */
349     if (rop->preop.clen % rop->imx.ncomp) {
350     fprintf(stderr, "%s: -c must have N x %d coefficients\n",
351     rop->inspec, rop->imx.ncomp);
352     return(0);
353     }
354     if (rop->preop.nsf > 0) { /* scale transform, instead */
355     if (rop->preop.nsf == 1) {
356     for (i = rop->preop.clen; i--; )
357     rop->preop.cmat[i] *= rop->preop.sca[0];
358     } else if (rop->preop.nsf*rop->imx.ncomp != rop->preop.clen) {
359     fprintf(stderr, "%s: -s must have one or %d factors\n",
360     rop->inspec,
361     rop->preop.clen/rop->imx.ncomp);
362     return(0);
363     } else {
364     for (i = rop->preop.nsf; i--; )
365     for (j = rop->imx.ncomp; j--; )
366     rop->preop.cmat[i*rop->imx.ncomp+j]
367     *= rop->preop.sca[i];
368     }
369     }
370     rop->preop.nsf = 0; /* now folded in */
371     if (!split_input(rop)) /* get our own struct */
372     return(0);
373     rop->rmp->ncomp = rop->preop.clen / rop->imx.ncomp;
374     if ((rop->rmp->ncomp > 3) & (rop->rmp->dtype <= DTspec)) {
375     rop->rmp->dtype = DTfloat; /* probably not actual spectrum */
376     memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
377     }
378     } else if (rop->preop.nsf > 0) { /* else use scalar(s)? */
379     if (rop->preop.nsf == 1) {
380     for (i = rop->rmp->ncomp; --i; )
381     rop->preop.sca[i] = rop->preop.sca[0];
382     rop->preop.nsf = rop->rmp->ncomp;
383     } else if (rop->preop.nsf != rop->rmp->ncomp) {
384     fprintf(stderr, "%s: -s must have one or %d factors\n",
385     rop->inspec, rop->rmp->ncomp);
386     return(0);
387     }
388     }
389     return(1);
390     }
391    
392     static int
393     apply_op(RMATRIX *dst, const RMATRIX *src, const RUNARYOP *ro)
394     {
395     if (ro->clen > 0) {
396     RMATRIX *res = rmx_transform(src, dst->ncomp, ro->cmat);
397     if (!res) {
398     fputs("Error in call to rmx_transform()\n", stderr);
399     return(0);
400     }
401     if (!rmx_transfer_data(dst, res, 0))
402     return(0);
403     rmx_free(res);
404     } else if (dst != src)
405 greg 2.9 memcpy(dst->mtx, src->mtx, rmx_array_size(dst));
406 greg 2.1 if (ro->nsf == dst->ncomp)
407     rmx_scale(dst, ro->sca);
408     return(1);
409     }
410    
411     static int
412     open_input(ROPMAT *rop)
413     {
414     int outtype;
415    
416     if (!rop || !rop->inspec || !rop->inspec[0])
417     return(0);
418     if (rop->inspec == stdin_name)
419     rop->infp = stdin;
420     else if (rop->inspec[0] == '!')
421     rop->infp = popen(rop->inspec+1, "r");
422     else
423     rop->infp = fopen(rop->inspec, "rb");
424    
425     if (!rmx_load_header(&rop->imx, rop->infp)) {
426     fprintf(stderr, "Bad header from: %s\n", rop->inspec);
427     return(0);
428     }
429     return(get_component_xfm(rop));
430     }
431    
432     /* Return nominal wavelength associated with input component (return nm) */
433     static double
434     l_wavelength(char *nam)
435     {
436     double comp = argument(1);
437    
438     if ((comp < -.5) | (comp >= in_ncomp+.5)) {
439     errno = EDOM;
440     return(.0);
441     }
442     if (comp < .5) /* asking for #components? */
443     return(in_ncomp);
444    
445     if (in_ncomp == 3) { /* special case for RGB */
446     const int w0 = (int)(comp - .5);
447     return(mop[0].rmp->wlpart[w0] +
448     (comp-.5)*(mop[0].rmp->wlpart[w0+1] -
449     mop[0].rmp->wlpart[w0]));
450     }
451     return(mop[0].rmp->wlpart[0] + /* general case, even div. */
452     (comp-.5)/(double)in_ncomp *
453     (mop[0].rmp->wlpart[3] - mop[0].rmp->wlpart[0]));
454     }
455    
456     /* Return ith input with optional channel selector */
457     static double
458     l_chanin(char *nam)
459     {
460     double inp = argument(1);
461     int mi, chan;
462    
463     if ((mi = (int)(inp-.5)) < 0 || mi >= nmats) {
464     errno = EDOM;
465     return(.0);
466     }
467     if (inp < .5) /* asking for #inputs? */
468     return(nmats);
469    
470     if (nargum() >= 2) {
471     double cval = argument(2);
472     if (cval < .5 || (chan = (int)(cval-.5)) >= in_ncomp) {
473     errno = EDOM;
474     return(.0);
475     }
476     } else
477     chan = cur_chan;
478    
479     return(mop[mi].rmp->mtx[cur_col*in_ncomp + chan]);
480     }
481    
482     static int
483     initialize(RMATRIX *imp)
484     {
485     int i;
486     /* XXX struct is zeroed coming in */
487     setcolor(imp->cexp, 1.f, 1.f, 1.f);
488     for (i = 0; i < nmats; i++) { /* open each input */
489     int restype;
490     if (!open_input(&mop[i]))
491     return(0);
492     restype = mop[i].rmp->dtype;
493     if (!imp->dtype || (restype = rmx_newtype(restype, imp->dtype)) > 0)
494     imp->dtype = restype;
495     else
496     fprintf(stderr, "%s: warning - data type mismatch\n",
497     mop[i].inspec);
498     if (!i) {
499     imp->ncols = mop[0].rmp->ncols;
500     imp->ncomp = mop[0].rmp->ncomp;
501     memcpy(imp->wlpart, mop[0].rmp->wlpart, sizeof(imp->wlpart));
502 greg 2.2 } else if ((mop[i].rmp->ncols != imp->ncols) |
503     (mop[i].rmp->ncomp != imp->ncomp) |
504     ((in_nrows > 0) & (mop[i].rmp->nrows > 0) &
505     (mop[i].rmp->nrows != in_nrows))) {
506 greg 2.1 fprintf(stderr, "%s: mismatch in size or #components\n",
507     mop[i].inspec);
508     return(0);
509     } /* XXX should check wlpart? */
510 greg 2.2 if (in_nrows <= 0)
511     in_nrows = imp->nrows = mop[i].rmp->nrows;
512 greg 2.1 } /* set up .cal environment */
513     esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
514     esupport &= ~(E_OUTCHAN|E_INCHAN);
515     varset("PI", ':', M_PI);
516     varset("nfiles", ':', nmats);
517     varset("nrows", ':', in_nrows);
518     varset("ncols", ':', in_ncols);
519     varset("ncomp", ':', in_ncomp);
520     varset("R", ':', 1.);
521     varset("G", ':', 2.);
522     varset("B", ':', 3.);
523     funset("wl", 1, ':', l_wavelength);
524     funset("ci", 1, '=', l_chanin);
525     scompile("ri(i)=ci(i,R);gi(i)=ci(i,G);bi(i)=ci(i,B)", NULL, 0);
526     return(1);
527     }
528    
529     static void
530     output_headinfo(FILE *fp)
531     {
532     int i;
533    
534     for (i = 0; i < nmats; i++) {
535     const char *cp = mop[i].imx.info;
536     fputs(mop[i].inspec, fp);
537     fputs(":\n", fp);
538     if (!cp) continue;
539     while (*cp) {
540     if (*cp == '\n') {
541     cp++; /* avoid inadvertant terminus */
542     continue;
543     }
544     fputc('\t', fp); /* indent this input's info */
545     do
546     putc(*cp, fp);
547     while (*cp++ != '\n');
548     }
549     }
550     }
551    
552     static int
553 greg 2.10 spawned_children(int np)
554 greg 2.1 {
555 greg 2.10 int i, rv;
556    
557     #if defined(_WIN32) || defined(_WIN64)
558     if (np > 1) {
559     fputs("Warning: only one process under Windows\n", stderr);
560     np = 1;
561     } else
562     #endif
563     if ((in_nrows > 0) & (np > in_nrows))
564     np = in_nrows;
565     /* we'll be doing a row at a time */
566 greg 2.1 for (i = 0; i < nmats; i++) {
567 greg 2.10 mop[i].imx.nrows = 1;
568 greg 2.1 if (!rmx_prepare(&mop[i].imx))
569     goto memerror;
570     if (mop[i].rmp != &mop[i].imx) {
571     mop[i].rmp->nrows = 1;
572     if (!rmx_prepare(mop[i].rmp))
573     goto memerror;
574     }
575     }
576 greg 2.10 /* prep output row buffer */
577     if (mcat || mop[nmats].preop.clen > 0) {
578     if (!split_input(&mop[nmats])) /* need separate buffer */
579     return(0);
580     if (mop[nmats].preop.clen > 0)
581     mop[nmats].rmp->ncomp = mop[nmats].preop.clen /
582     mop[nmats].imx.ncomp;
583     mop[nmats].rmp->nrows = 1;
584     if (!mcat | !mcat_last && !rmx_prepare(mop[nmats].rmp))
585 greg 2.1 goto memerror;
586     }
587 greg 2.10 mop[nmats].imx.nrows = 1;
588     if (!rmx_prepare(&mop[nmats].imx))
589     goto memerror;
590     if (np <= 1) { /* single process return point */
591     #ifdef getc_unlocked
592     for (i = 0; i < nmats; i++)
593     flockfile(mop[i].infp);
594     flockfile(stdout);
595     #endif
596     return(0);
597     }
598 greg 2.13 pgid = setpgrp(); /* set process group ID */
599 greg 2.12 signal(SIGIO, on_sigio);
600 greg 2.10 fflush(stdout); /* flush header & spawn children */
601     cproc = (SUBPROC *)malloc(sizeof(SUBPROC)*np);
602     if (!cproc)
603     goto memerror;
604     nchildren = np;
605     for (i = 0; i < np; i++) {
606     cproc[i].flags = PF_FILT_OUT;
607     cproc[i].w = dup(1);
608     cproc[i].r = 0;
609     cproc[i].pid = -1;
610     rv = open_process(&cproc[i], NULL);
611     if (rv <= 0) break;
612     }
613     if (rv > 0)
614     return(1); /* parent return value */
615     if (rv < 0) {
616     perror("fork");
617     exit(1);
618     }
619 greg 2.11 inchild = i; /* our child index */
620     while (i-- > 0) /* don't share siblings' pipes */
621     close(cproc[i].w);
622 greg 2.10 fpurge(stdin); /* discard previous matrix input */
623     #ifdef getc_unlocked
624     flockfile(stdin);
625     #endif
626     for (i = 0; i < nmats; i++) {
627     if (mop[i].infp != stdin)
628 greg 2.11 fclose(mop[i].infp); /* ! pclose() */
629 greg 2.10 mop[i].infp = stdin;
630     mop[i].imx.dtype = DTdouble;
631     }
632     return(0); /* child return */
633     memerror:
634     fputs("Out of memory in spawned_children()\n", stderr);
635     exit(1);
636     }
637    
638     static int
639     parent_loop()
640     {
641     FILE **outfp = (FILE **)malloc(nchildren*sizeof(FILE *));
642     int i;
643    
644     if (!outfp) goto memerror;
645     for (i = 0; i < nchildren; i++) {
646     outfp[i] = fdopen(cproc[i].w, "w");
647     if (!outfp[i]) goto memerror;
648     #ifdef getc_unlocked
649     flockfile(outfp[i]);
650     #endif
651     }
652     #ifdef getc_unlocked
653     for (i = 0; i < nmats; i++)
654     flockfile(mop[i].infp);
655     #endif
656     for (cur_row = 0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row++) {
657     FILE *ofp = outfp[cur_row % nchildren];
658 greg 2.11 for (i = 0; i < nmats; i++)
659 greg 2.10 if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
660     if (cur_row > in_nrows) /* unknown #input rows? */
661     break;
662     fprintf(stderr, "%s: read error at row %d\n",
663     mop[i].inspec, cur_row);
664     return(0);
665     }
666 greg 2.11 if (i < nmats)
667     break;
668     for (i = 0; i < nmats; i++)
669 greg 2.10 if (!rmx_write_data(mop[i].imx.mtx, mop[i].imx.ncomp,
670 greg 2.12 mop[i].imx.ncols, DTdouble, ofp))
671 greg 2.10 return(0);
672     if (fflush(ofp) == EOF)
673     return(0);
674     }
675     for (i = 0; i < nchildren; i++)
676     fclose(outfp[i]);
677     free(outfp);
678     i = close_processes(cproc, nchildren);
679     free(cproc); cproc = NULL;
680     if (i < 0) {
681     fputs("Warning: missing child in parent_loop()\n", stderr);
682     return(1);
683     }
684     if (i > 0) {
685     fprintf(stderr, "Child exited with status %d\n", i);
686     return(0);
687     }
688     return(1);
689     memerror:
690     fputs("Out of memory in parent_loop()\n", stderr);
691     exit(1);
692     }
693    
694     static int
695     combine_input()
696     {
697     const int row0 = (inchild >= 0)*inchild;
698     const int rstep = nchildren + !nchildren;
699     ROPMAT *res = &mop[nmats];
700     int set_r, set_c;
701     RMATRIX *tmp = NULL;
702     int co_set;
703 greg 2.13 sigset_t iomask;
704 greg 2.10 int i;
705    
706 greg 2.1 if (mcat && mcat_last &&
707     !(tmp = rmx_alloc(1, res->imx.ncols, res->rmp->ncomp)))
708     goto memerror;
709     /* figure out what the user set */
710     co_set = fundefined("co");
711     if (!co_set)
712     co_set = -vardefined("co");
713     if (!co_set & (in_ncomp == 3) && vardefined("ro") &&
714     vardefined("go") && vardefined("bo")) {
715     scompile("co(p)=select(p,ro,go,bo)", NULL, 0);
716     co_set = 1;
717     }
718     if (co_set) { /* set if user wants, didn't set */
719     set_r = varlookup("r") != NULL && !vardefined("r");
720     set_c = varlookup("c") != NULL && !vardefined("c");
721     } else /* save a little time */
722     set_r = set_c = 0;
723 greg 2.13
724     sigemptyset(&iomask); /* read/process row-by-row */
725     sigaddset(&iomask, SIGIO);
726 greg 2.10 for (cur_row = row0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row += rstep) {
727 greg 2.1 RMATRIX *mres = NULL;
728 greg 2.14 if (inchild >= 0) sigprocmask(SIG_BLOCK, &iomask, NULL);
729     for (i = 0; i < nmats; i++)
730 greg 2.1 if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
731 greg 2.5 if (cur_row > in_nrows) /* unknown #input rows? */
732 greg 2.10 break;
733 greg 2.1 fprintf(stderr, "%s: read error at row %d\n",
734     mop[i].inspec, cur_row);
735     return(0);
736     }
737 greg 2.14 if (inchild >= 0) sigprocmask(SIG_UNBLOCK, &iomask, NULL);
738 greg 2.10 if (i < nmats)
739     break;
740     for (i = 0; i < nmats; i++)
741 greg 2.1 if (!apply_op(mop[i].rmp, &mop[i].imx, &mop[i].preop))
742     return(0);
743     if (set_r) varset("r", '=', cur_row);
744     for (cur_col = 0; cur_col < in_ncols; cur_col++) {
745     if (set_c) varset("c", '=', cur_col);
746     for (cur_chan = 0; cur_chan < in_ncomp; cur_chan++) {
747     const int ndx = cur_col*in_ncomp + cur_chan;
748     eclock++;
749     if (!co_set) { /* just summing elements? */
750     res->imx.mtx[ndx] = 0;
751     for (i = nmats; i--; )
752     res->imx.mtx[ndx] += mop[i].rmp->mtx[ndx];
753     } else if (co_set > 0) {
754     double dchan = cur_chan+1;
755     res->imx.mtx[ndx] = funvalue("co", 1, &dchan);
756     } else
757     res->imx.mtx[ndx] = varvalue("co");
758     }
759     } /* final conversions */
760     if (!mcat) {
761     if (!apply_op(res->rmp, &res->imx, &res->preop))
762     return(0);
763     } else if (mcat_last) {
764     if (!apply_op(tmp, &res->imx, &res->preop))
765     return(0);
766     mres = rmx_multiply(tmp, mcat);
767     if (!mres)
768     goto multerror;
769     if (!rmx_transfer_data(res->rmp, mres, 0))
770     return(0);
771     } else /* mcat && !mcat_last */ {
772     mres = rmx_multiply(&res->imx, mcat);
773     if (!mres)
774     goto multerror;
775     if (!apply_op(res->rmp, mres, &res->preop))
776     return(0);
777     }
778     rmx_free(mres); mres = NULL;
779 greg 2.13 if (inchild >= 0) { /* children share stdout */
780 greg 2.12 while (nr_out < cur_row)
781     pause(); /* wait for our turn */
782 greg 2.13 sigprocmask(SIG_BLOCK, &iomask, NULL);
783     }
784 greg 2.1 if (!rmx_write_data(res->rmp->mtx, res->rmp->ncomp,
785 greg 2.10 res->rmp->ncols, res->rmp->dtype, stdout))
786     return(0);
787 greg 2.12 if (inchild >= 0) { /* flush and notify group */
788     if (fflush(stdout) == EOF)
789     return(0);
790 greg 2.13 sigprocmask(SIG_UNBLOCK, &iomask, NULL);
791     killpg(pgid, SIGIO); /* increments everyone's nr_out */
792 greg 2.12 }
793 greg 2.1 }
794 greg 2.10 return(inchild >= 0 || fflush(stdout) != EOF);
795 greg 2.1 memerror:
796     fputs("Out of buffer space in combine_input()\n", stderr);
797     return(0);
798     multerror:
799     fputs("Unexpected matrix multiply error in combine_input()\n", stderr);
800     return(0);
801     }
802    
803     static int
804     get_factors(double da[], int n, char *av[])
805     {
806     int ac;
807    
808     for (ac = 0; ac < n && isflt(av[ac]); ac++)
809     da[ac] = atof(av[ac]);
810     return(ac);
811     }
812    
813     static void
814     resize_inparr(int n2alloc)
815     {
816     int i;
817    
818 greg 2.6 if (n2alloc == nall)
819     return;
820     for (i = nall; i > n2alloc; i--) {
821 greg 2.1 rmx_reset(&mop[i].imx);
822     if (mop[i].rmp != &mop[i].imx)
823     rmx_free(mop[i].rmp);
824     }
825     mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
826     if (mop == NULL) {
827     fputs("Out of memory in resize_inparr()\n", stderr);
828     exit(1);
829     }
830 greg 2.6 if (n2alloc > nall)
831     memset(mop+nall, 0, (n2alloc-nall)*sizeof(ROPMAT));
832 greg 2.1 nall = n2alloc;
833     }
834    
835     /* Load one or more matrices and operate on them, sending results to stdout */
836     int
837     main(int argc, char *argv[])
838     {
839    
840     int outfmt = DTfromHeader;
841     const char *defCsym = NULL;
842     int echoheader = 1;
843     int stdin_used = 0;
844 greg 2.10 int nproc = 1;
845 greg 2.1 const char *mcat_spec = NULL;
846     int n2comp = 0;
847     uby8 comp_ndx[128];
848     int i;
849     /* get starting input array */
850     mop = (ROPMAT *)calloc(nall=2, sizeof(ROPMAT));
851     /* get options and arguments */
852     for (i = 1; i < argc; i++)
853     if (argv[i][0] != '-' || !argv[i][1]) {
854     if (argv[i][0] == '-') {
855     if (stdin_used++) goto stdin_error;
856     mop[nmats].inspec = stdin_name;
857     } else
858     mop[nmats].inspec = argv[i];
859     if (!mop[nmats].preop.csym)
860     mop[nmats].preop.csym = defCsym;
861     if (++nmats >= nall)
862     resize_inparr(nmats + (nmats>>2) + 2);
863     } else {
864     int n = argc-1 - i;
865     switch (argv[i][1]) { /* get option */
866     case 'w':
867     nowarn = !nowarn;
868     break;
869     case 'h':
870     echoheader = !echoheader;
871     break;
872 greg 2.10 case 'n':
873     nproc = atoi(argv[++i]);
874     if (nproc <= 0)
875     goto userr;
876     break;
877 greg 2.1 case 'e':
878     if (!n) goto userr;
879     comp_ndx[n2comp++] = i++;
880     break;
881     case 'f':
882     switch (argv[i][2]) {
883     case '\0':
884     if (!n) goto userr;
885     comp_ndx[n2comp++] = i++;
886     break;
887     case 'd':
888     outfmt = DTdouble;
889     break;
890     case 'f':
891     outfmt = DTfloat;
892     break;
893     case 'a':
894     outfmt = DTascii;
895     break;
896     case 'c':
897     outfmt = DTrgbe;
898     break;
899     default:
900     goto userr;
901     }
902     break;
903     case 's':
904     if (n > MAXCOMP) n = MAXCOMP;
905     i += mop[nmats].preop.nsf =
906     get_factors(mop[nmats].preop.sca,
907     n, argv+i+1);
908     if (mop[nmats].preop.nsf <= 0) {
909     fprintf(stderr, "%s: -s missing arguments\n",
910     argv[0]);
911     goto userr;
912     }
913     break;
914     case 'C':
915     if (!n || isflt(argv[i+1]))
916     goto userr;
917     defCsym = mop[nmats].preop.csym = argv[++i];
918     mop[nmats].preop.clen = 0;
919     mcat_last = 0;
920     break;
921     case 'c':
922     if (n && !isflt(argv[i+1])) {
923     mop[nmats].preop.csym = argv[++i];
924     mop[nmats].preop.clen = 0;
925 greg 2.15 mcat_last = 0;
926 greg 2.1 break;
927     }
928     if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
929     i += mop[nmats].preop.clen =
930     get_factors(mop[nmats].preop.cmat,
931     n, argv+i+1);
932     if (mop[nmats].preop.clen <= 0) {
933     fprintf(stderr, "%s: -c missing arguments\n",
934     argv[0]);
935     goto userr;
936     }
937     mop[nmats].preop.csym = NULL;
938     mcat_last = 0;
939     break;
940     case 'm':
941     if (!n) goto userr;
942     if (argv[++i][0] == '-' && !argv[i][1]) {
943     if (stdin_used++) goto stdin_error;
944     mcat_spec = stdin_name;
945     } else
946     mcat_spec = argv[i];
947     mcat_last = 1;
948     break;
949     default:
950     fprintf(stderr, "%s: unknown option '%s'\n",
951     argv[0], argv[i]);
952     goto userr;
953     }
954     }
955     if (!nmats) {
956     fprintf(stderr, "%s: need at least one input matrix\n", argv[0]);
957     goto userr;
958     }
959     resize_inparr(nmats+1); /* extra matrix at end for result */
960     mop[nmats].inspec = "trailing_ops";
961     /* load final concatenation matrix */
962     if (mcat_spec && !(mcat = rmx_load(mcat_spec, RMPnone))) {
963     fprintf(stderr, "%s: error loading concatenation matrix: %s\n",
964     argv[0], mcat_spec);
965     return(1);
966     }
967     /* get/check inputs, set constants */
968     if (!initialize(&mop[nmats].imx))
969     return(1);
970    
971     for (i = 0; i < n2comp; i++) /* user .cal files and expressions */
972     if (argv[comp_ndx[i]][1] == 'f') {
973     char *fpath = getpath(argv[comp_ndx[i]+1],
974     getrlibpath(), 0);
975     if (fpath == NULL) {
976     fprintf(stderr, "%s: cannot find file '%s'\n",
977     argv[0], argv[comp_ndx[i]+1]);
978     return(1);
979     }
980     fcompile(fpath);
981     } else /* (argv[comp_ndx[i]][1] == 'e') */
982     scompile(argv[comp_ndx[i]+1], NULL, 0);
983    
984     /* get trailing color transform */
985     if (!get_component_xfm(&mop[nmats]))
986     return(1);
987     /* adjust output dimensions and #components */
988     if (mcat) {
989     if (mop[nmats].imx.ncols != mcat->nrows) {
990     fprintf(stderr,
991     "%s: number of input columns does not match number of rows in '%s'\n",
992     argv[0], mcat_spec);
993     return(1);
994     }
995     if (mcat->ncomp != (mcat_last ? mop[nmats].rmp->ncomp : mop[nmats].imx.ncomp)) {
996     fprintf(stderr,
997     "%s: number of components does not match those in '%s'\n",
998     argv[0], mcat_spec);
999     return(1);
1000     }
1001     if (!split_input(&mop[nmats]))
1002     return(1);
1003     mop[nmats].rmp->ncols = mcat->ncols;
1004     }
1005     newheader("RADIANCE", stdout); /* write output header */
1006     if (echoheader)
1007     output_headinfo(stdout);
1008     printargs(argc, argv, stdout);
1009     fputnow(stdout);
1010     mop[nmats].rmp->dtype = rmx_write_header(mop[nmats].rmp, outfmt, stdout);
1011     if (!mop[nmats].rmp->dtype) {
1012     fprintf(stderr, "%s: unsupported output format\n", argv[0]);
1013     return(1);
1014     }
1015 greg 2.7 doptimize(1); /* optimize definitions */
1016 greg 2.11 if (spawned_children(nproc)) /* running in parent process? */
1017 greg 2.10 return(parent_loop() ? 0 : 1);
1018 greg 2.1 /* process & write rows */
1019 greg 2.10 return(combine_input() ? 0 : 1);
1020 greg 2.1 stdin_error:
1021     fprintf(stderr, "%s: %s used for more than one input\n",
1022     argv[0], stdin_name);
1023     return(1);
1024     userr:
1025     fprintf(stderr,
1026 greg 2.10 "Usage: %s [-h][-f{adfc}][-n nproc][-e expr][-f file][-s sf .. | -c ce ..] m1 .. -m mcat > mres\n",
1027 greg 2.1 argv[0]);
1028     return(1);
1029     }