ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/rcomb.c
Revision: 2.35
Committed: Thu Oct 30 16:46:49 2025 UTC (9 days, 16 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.34: +7 -2 lines
Log Message:
feat(rcomb): Added proper error message when input file cannot be opened

File Contents

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