ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.12
Committed: Tue May 21 17:39:17 2024 UTC (11 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.11: +22 -8 lines
Log Message:
fix(rcomb): Switched to group-signaling to make multi-processing more reliable

File Contents

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