ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.18
Committed: Thu May 23 19:29:41 2024 UTC (10 days, 8 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.17: +18 -6 lines
Log Message:
fix(rcomb): Avoid inappropriate byte-swapping and improved final fclose() wait

File Contents

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